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ABSTRACT 

The  main  goal  of  this  research  is  to  model  a  flexible  missile  with  structural 
flexibility  utilizing  the  Equivalent  Rigid  Link  System  (ERLS)  with  an  enhanced 
natural  mode  discretization.  Dynamic  analysis  of  the  flexible  missile  in 
2-Dimension  motion  is  presented. 

Computer  simulation  is  performed  where  the  pitch  angle  of  the  missile  is 
controlled  with  a  rigid-body  controller.  The  effects  of  increasing  payloads  and 
speed  to  the  system  performance  are  discussed. 
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THESIS  DISCLAIMER 

The  reader  is  cautioned  that  computer  programs  developed  in  this  research 
may  not  have  been  exercised  for  all  cases  of  interest.  While  every  effort  has  been 
made,  within  the  time  available,  to  ensure  that  the  programs  are  free  of  compu- 
tational and  logic  errors,  they  cannot  be  considered  validated.  Any  appUcation 
of  these  programs  without  additional  verification  is  at  the  risk  of  the  user. 


TABLE  OF  CONTENTS 

I.     INTRODUCTION 1 

A.  BACKGROUND    1 

B.  LITERATURE  REVIEW    1 

C.  OUTLINE    3 

IL     MODEL  FORMULATION  FOR  A  FLEXIBLE  MISSILE     4 

A.  KINEMATICS    .'. 4 

B.  KINETICS    8 

III.  THE  DEVELOPMENT  OF  A  RIGID-BODY  CONTROLLER 12 

IV.  COMPUTER  SIMULATION     16 

A.  SIMULATION  OBJECTIVES     16 

B.  SOLUTION  TECHNIQUE      16 

C.  THE  COMPUTER  SIMULATION  CODE   17 

D.  SIMULATION  RESULTS    17 

V.  SUMMARY    49 

A.  CONCLUSIONS    49 

B.  RECOMMENDATIONS     49 

APPENDIX  A.     DERIVATION  OF  THE  MODE  SHAPE  RESPONSE 
MATRIX  FOR  THE  FLEXIBLE  MISSILE    50 

APPENDIX  B.     DERIVATION  OF  THE  EQUATIONS  OF  MOTION 
FOR  THE  FLEXIBLE  MISSILE 57 


APPENDIX  C.     DERIVATION  OF  THE  GENERALIZED  FORCES  OF 
THE  EQUATIONS  OF  THE  MOTION    65 

APPENDIX  D.     COMPUTER  SIMULATION  CODE    67 

LIST  OF  REFERENCES     91 

V'' 

INITIAL  DISTRIBUTION  LIST     92 


VI 


LIST  OF  FIGURES 

Figure     1.  The  general  configuration  of  the  missile    4 

Figure     2.  Equivalent  Rigid  Link  System  (ERLS)     5 

Figure     3.  Applied  Generalized  Forces  on  Flexible  Missile     9 

Figure    4.  The  flexible  missile  with  rigid-body  controller 12 

Figure     5.  The  desired  and  actual  trajectory  for  the  rigid  missile  with  rigid- 
body  controller    19 

Figure     6.  The  control  angle  for  the  rigid  missile  with  rigid-body  controller     20 
Figure     7.  The  large  motion  position  X  for  the  rigid  missile  with  rigid-body 

controller     21 

Figure     8.  The  large  motion  position  Y  for  the  rigid  missile  with  rigid-body 

controller     22 

Figure     9.  The  desired  and  actual  trajectory  for  the  flexible  missile  with 

rigid-body  controller    23 

Figure  10.  The  control  angle  for  the  flexible  missile  with  rigid-body  controller  24 
Figure  1 1.  The  small  motion  position  v(0)  for  the  flexible  missile  with  rigid- 
body  controller    25 

Figure  12.  The  small  motion  position  phi(O)  for  the  flexible  missile  with 

rigid-body  controller    26 

Figure  13.  The  large  motion  position  X  for  the  flexible  missile  with  rigid-body 

controller    27 

Figure  14.  The  large  motion  position  Y  for  the  flexible  missile  with  rigid-body 

controller     28 

Figure  15.  The  desired  and  actual  trajectory  for  the  flexible  missile  with 

rigid-body  controller  without  saturation  on  the  control  angle     .  .    30 
Figure  16.  The  control  angle  for  the  flexible  missile  with  rigid-body  controller 

without  saturation  on  the  control  angle     31 

Figure  17.  The  small  motion  position  v(0)  for  the  flexible  missile  rigid-body 

controller  without  saturation  on  the  control  angle 32 


vu 


Figure  18.  The  small  motion  position  phi(O)  for  the  flexible  missile  rigid-body 

controller  without  saturation  on  the  control  angle 33 

Figure  19.  The  large  motion  position  X  for  the  flexible  missile  rigid-body 

controller  without  saturation  on  the  control  angle 34 

Figure  20.  The  large  motion  position  Y  for  the  flexible  missile  rigid-body 

controller  without  saturation  on  the  control  angle 35 

Figure  21.  The  desired  and  actual  trajectory  for  the  flexible  missile  using 

rigid-body  controller  with  saturation  on  the  control  angle    37 

Figure  22.  The  control  angle  for  the  flexible  missile  using  rigid-body  control- 
ler with  saturation  on  the  control  angle     38 

Figure  23.  The  small  motion  position  v(0)  for  the  flexible  missile  using  rigid- 
body  controller  with  saturation  on  the  control  angle    39 

Figure  24.  The  small  motion  position  phi(O)  for  the  flexible  missile  using 

rigid-body  controller  with  saturation  on  the  control  angle    40 

Figure  25.  The  large  motion  position  X  for  the  flexible  missile  using  rigid- 
body  controller  with  saturation  on  the  control  angle    41 

Figure  26.  The  large  motion  position  Y  for  the  flexible  missile  using     42 

Figure  27.  The  desired  and  actual  trajectory  for  the  flexible  missile  using 

rigid-body  controller  with  increased  speed     43 

Figure  28.  The  control  angle  for  the  flexible  missile  using  rigid-body  control- 
ler with  increased  speed   44 

Figure  29.  The  small  motion  position  v(0)  for  the  flexible  missile  using  rigid- 
body  controller  with  increased  speed 45 

Figure  30.  The  small  motion  position  phi(O)  for  the  flexible  missile  using 

rigid-body  controller  with  increased  speed     46 

Figure  31.  The  large  motion  position  X  for  the  flexible  missile  using  rigid- 
body  controller  with  increased  speed     47 

Figure  32.  The  large  motion  position  Y  for  the  flexible  missile  using  rigid- 
body  controller  with  increased  speed     48 

Figure  33.  The  Bar  in  Flexure    50 

Figure  34.  Free  Body  Diagram  Corresponding  to  a  Bar  Element   51 


vm 


ACKNOWLEDGEMENTS 

I  would  like  to  extend  my  sincerest  thanks  and  appreciation  to  my  advisor, 
Professor  Liang-Wey  Chang  for  his  steadfast  support,  guidance  and  always  open 
door  policy.  I  would  also  like  to  thank  Professor  Harold  A.Titus  and  Laboratory 
Technician  Colin  Cooper  for  their  technical  support. 

I  would  especially  Hke  to  recognize  the  accomplishments  of  my  wife,  Ulfet. 
Her  steadfast  devotion,  encouragement,  support  and  confidence,  permitted  my 
absence  during  hours  of  derivations  and  computer  work  while  she  carried  the 
burden  of  household  and  family  work. 


IX 


I.     INTRODUCTION 

A.  BACKGROUND 

With  the  introduction  of  long  slender  missiles  such  as  the  Vanguard,  the 
Redstone,  and  various  ballistic  missiles,  the  problem  of  the  structural  flexibility 
became  severe  [Ref.  1].  Due  to  the  limited  thrust  available  from  rocket  engines, 
these  missiles  had  to  be  as  light  as  possible. This  meant  a  sacrifice  in  structural 
rigidity. 

Missile  flexure  causes  additional  aerodynamic  loads  which  in  turn  cause  ad- 
ditional flexure.  Coupling  occurs  between  the  elastic  modes  and  the  control  sys- 
tem as  the  control  system  gyros  sense  the  fiexure  motion  and  the  rigid  body 
motion.  It  has  become  necessary  to  actively  control  the  flexible  structure  and 
thereby  reduce  the  structural  loads  and  improve  the  \ehicle  response  such  as  po- 
sition, velocity  and  acceleration.  Reduced  structural  loads  will  also  offer  potential 
for  reduced  bending  stress  and  fatigue  problems. 

The  objective  of  this  thesis  is  to  develop  a  dynamic  model  for  a  flexible  missile 
and  study  the  dynamic  behavior  of  the  flexible  missile.  Simulation  is  a  valuable 
tool  in  the  design  of  new  missile  systems  and  in  the  modification  or  evaluation 
of  existing  systems.  A  missile  simulation  allows  the  engineer  to  evaluate  his  design 
without  the  expense  of  actually  building  and  flying  the  actual  missile.  System 
dynamics  can  be  investigated  through  simulation  with  a  substantial  savings  in 
time  and  expense  [Ref.  2]. 

B.  LITERATURE  REVIEW 

Flexible  missile  modeling  centers  on  the  relationship  between  the  large,  rigid- 
body  motion  and  the  small  motions  due  to  structural  flexibility. 

Jenkins  [Ref.  2]  expresses  some  techniques  used  in  deriving  the  equations  of 
motion  of  a  rigid  missile  for  a  six  degree-of-freedom  (6-DOF)  simulation.  The 
rigid  missiles  are  characterized  by  their  larger  size  and  low  ratio  of  payload  to 
total  weight.  Rigid  missile  dynamic  equations  were  developed  using  the 
Newton-Euler  Method.    The  moments  and  forces  along  with  the  mass  and  mo- 


ments  of  inertia  are  assumed  to  be  known  in  the  body  coordinate  frame.  The 
transformation  between  global  and  local  coordinate  frame  is  achieved  with  a 
non-homogenous  coordinate  transformation  matrix  [Ref.  3:  pp.  342].  The  result- 
ing rigid-body  equations  of  motion  produced  the  understanding  for  the  derivation 
of  the  dynamic  equations  of  the  flexible  missile. 

The  Equivalent  Rigid  Link.  System  [Ref.  4]  describes  the  large-motion 
kinematics  of  the  system  by  the  ERLS  motion  and  the  small  motion  relative  to 
the  ERLS.  The  application  of  the  finite  element  techniques  and  Lagrangian  dy- 
namics produces  two  sets  of  coupled,  nonhnear,  ordinary  differential  equations 
of  motion,  of  which  one  set  is  for  the  large  motions  and  the  other  set  for  the  small 
motions.  The  small  motion  is  described  by  the  superposition  of  vibration  modes. 
The  modes  of  the  vibration  of  the  flexible  bar  was  derived  with  the  simple-beam 
theory  [Ref.  5:  pp.221].  In  simple  beam  theory,  it  is  assumed  that  the  rotation  of 
the  element  is  insignificant  compared  to  the  vertical  translation  and  the  shear 
deformation  is  small  relative  to  the  bending  deformation.  This  assumption  is  valid 
if  the  ratio  between  the  length  of  the  bar  and  its  height  is  relatively  large  (more 
than  10).  The  set  of  large  motion  equations  are  non-hnear  in  both  the  large  and 
small  motion  variables  while  the  set  of  small  motion  equations  are  linear  in  the 
small  motion  variable  and  non-linear  in  the  large  motion  variables.  A  solution 
technique  called  the  Sequential  Integration  Method  [Ref.  6]  was  developed  which 
allows  efficient  simulations  of  systems  with  inertia  coupled  motions  having  non- 
linear slow  motion  (large  motion)  with  hnear  fast  motion  (small  motion).  The 
ERLS  model  presents  a  complete,  efficient  dynamic  model  able  to  describe  large 
motion,  small  motion  and  their  coupling. 

An  ERLS  model  of  a  flexible  spacecraft  boom  was  developed  and  a  computer 
simulation  was  performed  [Ref.  7].  The  equations  of  motion  were  solved  using 
the  Sequential  Integration  Method.  A  spatial  finite  discretization  of  the  boom 
structure  and  the  application  of  an  assumed  polynomial  modal  response  were 
utilize^  in  the  approximate  solution  to  the  equations  of  motion.  This  work  was 
the  basis  work  for  the  modeling  of  the  flexible  missiles. 


Ganon  [Ref.  8]  performed  an  experimental  validation  of  the  ERLS  dynamic 
model  in  a  vertical  plane  motion.  The  small  motion  was  modeled  using  a  shape 
matrix  derived  from  superposition  of  natural  modes.  The  agreement  between  the 
simulation  results  and  the  experimental  results  justify  the  application  of  the 
ERLS  using  a  natural-mode  shape  matrix  to  model  the  dynamic  response  of 
flexible  missile. 

C.     OUTLINE 

In  this  study,  the  ERLS  dynamic  model  is  used  to  derive  the  system  equations 
of  motion  for  a  flexible  missile  in  2-D  motion.  Dynamic  response  is  predicted  by 
solving  the  equations  of  motion  using  the  Sequential  Integration  Method.  The 
application  of  an  assumed  natural  mode  shape  and  the  spatial  finite  element 
discretization  of  the  missile  provide  an  approximate  solution.  Computer  simu- 
lation for  the  flexible  missile  is  performed  using  MATLAB  on  the  MACINTOSH 
computer.  A  rigid  body  controller  is  included  in  the  simulation  to  control  the 
pitch  angle  of  the  flexible  missile. 

The  ERLS  dynamic  model  of  the  flexible  missile  is  presented  in  Chapter  Two. 
A  rigid  body  controller  for  the  flexible  missile  is  described  in  Chapter  Three.  The 
computer  simulation  methodology  and  simulation  results  are  presented  in  Chap- 
ter Four.  The  conclusions  and  recommendations  for  future  work  are  presented 
in  Chapter  Five. 


II.     MODEL  FORMULATION  FOR  A  FLEXIBLE  MISSILE 

A.     KINEMATICS 

In  our  model,  2-D  inertial  reference  frame,  i.e.,  X  and  Y,  is  used  for  the  global 
motion  as  shown  in  Fig.l.  The  body-fixed  coordinate  frame,  i.e.,  x  and  y,  is  se- 
lected to  describe  the  missile  local  motion. 


A 


Y 


x=Local  X  axis 
y=Local  y  axis 
X=Global  X  axis 
Y=G1obal  Y  axis     • 
L=M1ss11e  Length 
D=Miss11e  Diameter 


Figure  I.      The  general  configuration  of  the  missile 


The  following  assumptions  were  made  for  the  fiexible  missile: 

(1)  Material  density  is  constant  throughout  the  body  and  the  steel  was  chosen 
for  modeling. 

(2)  A  uniform  circular  cross  section  is  assumed. 

The  geometric  configuration  parameters  and  material  properties  of  the  fiexi- 
ble missile  are  listed  as 


Diameter  =  0.1 2  meter 

Length  =  4.0  meter 

Material  Density  =  7861.05  kg/m^ 

Young's  Modulus  =  2.0  x  10"  pascal 

The  concept  of  the  ERLS  is  applied  to  model  the  kinematics  of  flexible  mis- 
sile. The  main  idea  of  the  Equivalent  Rigid  Link  System  (Fig. 2)  is  to  separate  the 
kinematics  of  the  flexible  body  into  large  and  small  motions. 


^^ 

t                  ^^^fl 

y 

^^^ 

\            ^^ 

w[^^^^                            ^ 

XJ^T 

^  ®  v(o) 

=  Missile   Base   Deflection 

ivTy 

0(0) 

=   Missile    Base   Slope 

/  oVf 

v(L) 

=   Missile   Tip    Deflection 

/  i 

0(L) 

=     Missile  Tip  Slope 

Y      K^ 

R 

=  Global   Position  Vector 

A  _♦    r^wif 

of   Missile   Base 

Ro/   / 

Ro 

=  Absolute  Position  Vector 

/  /-*- 

w 

of  the  base  of  ERLS 

/   /  R 

^^ 

/  / 

r 

=  Local  Position  Vector 

// 

d 

=  Deformation  Vector 

11 

. 

=    Equivalent    Rigid    Link 
=    Tfieoretlcal    Missile 

/ 

f 

^ 

Position 

X 

\ 

Figure  2.      Equivalent  Rigid  Link  System  (ERLS) 


The  large  motion  of  the  missile  will  be  represented  by  missile  base  position  in 
X  direction  Xq,  missile  base  position  in  Y  direction  Fq,  and  missile  pitch  angle  6. 
The  small  motions  resulting  from  flexible  motion  are  measured  relative  to  the  lo- 
cal coordinate  frame  x,  y.  v(0)  and  <t>{0)  are  the  nodal  displacement  and  slope 
of  the  missile  base  respectively.  The  absolute  (global)  position  of  a  point  on  the 
flexible  missile  is  obtained  from  coordinate  transformations.  The  global  position 
(R)  of  any  point  position  can  be  defined  using  ERLS  in  terms  of  a  local  position 
vector  (f),  a  deformation  vector  (d),  and  a  coordinate  transformation  matrix  (W), 
i.e., 


R  =  W{y  +  d) 


(2-1) 


The  transformation  matrix  (W)  between  the  large  motion  and  small  motion  co- 
ordinate is 


W  = 


1         0  0 

Xq   cos(0)    sin(^) 
Yq    sin(^)    cos((9) 


(2-2) 


and  the  local  rigid  body  position  vector  is 


r  — 


1 

X 

0 


(2-3) 


The  deformation  vector  (d)  is  expressed  in  terms  of  a  nodal  displacement 
vector  U  and  a  shape  function  </>  as 


d  =<t)U 
where 


(2-4) 


(/  =  [v(0)    0(0)]' 


(2-5) 


The  shape  function  0  was  derived  utilizing  a  natural-mode  superposition  and 
a  finite  element  concept.  The  Finite  Element  Method  (FEM)  was  utilized  to 
discretize  the  flexible  body  displacements  and  assigning  the  nodes.  Displacements 
are  for  each  point  along  the  missile,  a  function  of  location  and  time,  and  it  is 
necessary  to  discretize  the  deformations  to  obtain  an  ordinary  differential 
equation.  The  natural  mode  shape  function  of  a  beam  is  used  to  represent  the 
flexural  motion  of  the  flexible  missile.  Only  the  first  two  mode  shapes  are  used. 
The  fiexible  missile  is  modeled  as  a  continuous  Euler- Bernoulli  free-free  beam, 
neglecting  shear  deformation  and  rotary  inertia  effects.  The  nodal  points  are  the 
base  points  of  the  flexible  missile.  Appendix  A  shows  the  mode  shape  function 
derivation.   In  this  case,  0  is  found  as, 

"0   O' 
0=    0    0  (2-6) 


"0 

o' 

0 

0 

a 

b_ 

where  a  and  b  are  defined  as  following, 

a  =  /'i(Ci(  cos  p^x  +  cosh  P^x)  +  ( sin  P^x  +  sinh  p^x)) 

+  F^{C2{  cos  P2X  +  cosh  ^2-^)  +  ( sin  P2X  +  sinh  p2x))  (2  -  7) 

b  =  F2{Ci{  cos  pyX  +  cosh  pyx)  +  ( sin  PiX  +  sinh  Pix)) 

+  /^t(C2(  cos  P2X  +  cosh  P2X)  +  (  sin  P2X  +  sinh  P2X))  (2  -  8) 

Lagrangian  dynamics  is  used  to  acquire  the  equations  of  motion  and  the 
kinetic  energy  of  the  system  will  be  needed  in  the  development  of  the  Lagrangian 
expression.  The  absolute  velocity  is  listed  as  follows, 

R  =  W{T  +  d)-{-  Wd  (2-9) 


B.     KINETICS 

Lagrangian  Dynamics  is  a  systematic  way  to  derive  equations  of  motion  for 
complex  systems  like  flexible  missiles.   Lagrange's  equation  is  written  as, 

l(¥'-¥  +  ^  =  a         «=1,2,...,«)  (2-10) 

where 

K.E.  =  Kinetic  Energy 

P.E.  =  Potential  Energy 

q,  =  Generalized  coordinates 

O,  =  Generalized  forces 

n  =  Number  of  degrees  of  freedom 
The  generalized  coordinates  are  chosen  to  be, 

q  =  lXo   Yo  6  v(0)  0(0)]^  (2-11) 

The  kinetic  energy  of  the  system  is  defined  as  follows, 

A:.£.  =  yj^^^Jm  (2-  12) 

By  substituting  Eq.(2-1)  into  Eq.   (2-12),  the  kinetic  energy  is  written  as 


+  U'^(I)'^W'^IV(1)U  +  2U'^4)'^W'^W4)U  +  u'^4)'^w'^W4>U^)dm      ■  (2  -  13) 

The  potential  energy  of  the  system  includes  the  strain  energy  of  the  flexible 
missile  and  the  gravitational  potential  energy. 

P.E.  =  y  lu'^r^CrUdx  -  iR^gdm  (2-14) 


where 

r  =  Spatial  derivative  of  the  shape  function 

C  =  Rigidity  Matrix 

g  =  Gravitational  acceleration  vector 

Appendix  B  includes  the  development  of  the  equations  of  motions  using 
Lagrange  equations. 

Generalized  forces  will  be  found  by  virtual  work  principle.  It  is  assumed  that 
the  only  force  is  the  thrust  force  which  is  applied  to  the  base  of  the  flexible  missile 
as  shown  in  Fig. 3.  Other  forces  like  aerodynamic  and  damping  forces  are  neg- 
lected. 
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Figure  3.      Applied  Forces  on  Flexible  Missile 


Tlrt^hrust  force  is  composed  of  two  components  i.e., 
7^=  rcos(5  +  cD(0)) 


(2-15) 


T=  rsin(d  +  O(0)) 


(2-16) 


and  a  moment  applied  on  the  missile  is, 
M0  =  -rcos(5  +  (I>(O))v(O) 


(2-17) 


By  applying  virtual  work  principle,  the  generalized  forces  of  flexible  missile  is 
found  below. 


F>  = 


rcos(0)  -  T5  sin(0)  -  70(0)  sin(O) 

rsin(O)  +  n  cos(0)  +rO(0)  cos(0) 

-7v(0) 


(2-18) 


Fu  = 


Td  +  70(0) 
-Tv{0) 


(2-19) 


where 

T  =  Force 

d  =  Control  Angle  being  assumed  small 

F^  =  Large  motion  generalized  force  vector 

F^  =  Small  motion  generalized  force  vector 

Appendix  C  shows  the  derivation  of  the  generalized  forces  using  the  Virtual 
Work  Principle. 

The  derivation  of  the  equations  of  motion  from  the  Lagrangian  formulation 
yields  two  sets  of  coupled  equations.  One  set  describes  the  large  motions  and  the 
another  set  describes  the  small  motions.  These  two  sets  of  equations  are  non- 
linear, coupled,  second-order,  ordinary  differential  equations  of  the  form. 


M,,  ^   +  M^,  f/    =  F, 


(2  -  20) 
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M,^  ^  +  M,„  U  +  G„  U  +  K„  U    =  f,  (2  -  21) 

where: 

Al^q    =     3x3    effective  inertia  matrix  for  large  motions 

Mq„    =     3x2     coupled  inertia  matrix  of  the  small  motion  effect  on  large 
motions 

Fq  =     3x1     load  vector  for  the  large  motions 

M„^   =     2x3     coupled  inertia  matrix  of  the  large  motion  effect  on  small 
motions 

M„„    =     2x2    effective  inertia  matrix  for  small  motions 

G„  =     2x2    gyroscopic  matrix 

K„  =     2x2    stiffness  matrix 

F„  =     2x1     load  vector  for  the  small  motions 

^     =     3x1     vector,  generalized  coordinates  of  the  large  motions 

U    =     2x1     vector,  generalized  coordinates  of  the  small  motions 
The  detailed  development  of  the  equations  of  motion  and  definitions  of  the  terms 
are  described  in  Appendix  B. 
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III.     THE  DEVELOPMENT  OF  A  RIGID-BODY  CONTROLLER 

In  our  research,  the  controller  is  developed  based  on  the  rigid-body  assump- 
tion. The  purpose  of  including  the  rigid-body  controller  is  to  perform  computer 
simulation  and  study  the  dynamic  behavior  of  the  flexible  missile.  The  pitch  an- 
gle of  the  missile  is  controlled  by  the  rigid-body  controller  which  is  a  single  input 
single  output  (SISO)  controller.  A  general  configuration  of  the  flexible  missile 
and  rigid-body  controller  are  shown  in  Fig.4. 
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rigure  4.      The  flexible  missile  with  rigid-body  controller 

The  pitch  angle  is  the  control  variable  and  the  desired  states  are  desired  pitch 
angle  Bj,  desired  angular  velocity  6j,  and  desired  angular  acceleration  6^.  The 
control  angle  (5)  is  assumed  small. 

The  desired  error  function  is  defined  as 


e  +  Kyb  +  Kpt  =  0 


(3-1) 


where, 


12 


Kp  =  The  position  feedback  gain. 
K^  =  The  velocity  feedback  gain. 
e  =  6  -  dj  =  The  error  in  the  position. 
e  =  6  —  dj  =  The  error  in  velocity. 
t  =  d  —  dj  =  The  error  in  acceleration. 

The  control  design  begins  with  the  equation  of  rigid-body  (large  motion)  that 
is  obtained  by  modifying  and  expanding  Eq.(2-20)  as 


^n     ^n 


A/21      ^^ 


22 


r^i 

7r 

+ 

'h; 

+  6 

'b,' 

L'^2j 

UJ 

W\ 

Wi 

(3-2) 


Where, 


V/ii  = 


V/    (1,1)     .VU1,2) 


w 


^qq 


qq\ 


M,A2,\)     ALJ2,2) 


qq\ 


(3-3) 


Mn  = 


^V/,,(l,3) 
A/,,(2,3) 


(3-4) 


^h^  =  \.^^,,{^^)    AC(3,2)] 


qq 


qq\ 


^hi  =  [M,,(3,3)] 


qq\ 


(3-5) 
(3-6) 


'?!  = 


^0 


(3-7) 


ni  =  [^] 


(3-8) 


/  =  Gravity,  centrifugal  and  coriolis  forces. 
J2  =  Gravity,  centrifugal  and  coriolis  forces. 
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h,= 


Tcos{d) 
Tsm{e) 


(3-9) 


h2  =  [0] 


^,  = 


-Tsm{d) 
Tcos{d) 


(3  -  10) 


^2  =  CO] 


Eq.(3-2)  can  be  rewritten  in  a  tensor  form  as 


A/n^i  +  iV/i2^2  =/i  +  -^1  +  b^5 


(3-11) 


M2i^l4-M22?2=/2 


From  Eq.(3-ll), 


^1  =  A/f/ 


/i  +  /Zi  +  b^d  -  Mi2?2 


Substituting  Eq.(3-13)  into  Eq.(3-12),  we  find, 


(3-12) 


(3-13) 


Al]2  =  F  -  BS 


Where, 


-1 


A  =  M22  -  M2iMn  Mi2 


F=f2-M2iMu'[f,  +  h,] 


B  =  M2iMn% 


(3  -  14) 


(3-15) 
(3  -  16) 

(3-17) 
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Eq.(3-1)  can  be  rewritten  as 

(d  -  e,)  +  KXd  -  e^)  +  Kp{Q  -e^)  =  o  (3  -  i8) 

I 

By  substituting  Eq.(3-14)  into  Eq.(3-18),  we  find  the  control  angle  (8) 

s  =  -B-'[A(d^  -  KXe  -  e^)  -  Kp{e  -  e^))  -f]  (3  -  19) 

The  control  angle  {8)  is  a  function  of  the  pitch  angle  {6)  ,  desired  pitch  angle 
(dj),  angular  velocity  (6),  desired  angular  velocity  (dj),  desired  angular  acceler- 
ation (6j),  the  position  feedback  gain  Kp  ,  the  velocity  feedback  gain  K^,  and  the 
matrices  i.e.,  A,  B,  F.  Since  the  thrust  force  T  is  included  in  F  and  B,  the  mag- 
nitude of  the  thrust  force  will  thus  affect  the  control  angle  6  .  K^  and  Kp  are  ad- 
justable to  obtain  the  desired  response  of  the  pitch  angle  d. 
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IV.     COMPUTER  SIMULATION 

A.  SIMULATION  OBJECTIVES 

In  literature,  the  application  of  forces  has  been  limited  to  gravity,  aerodyna- 
mic forces  and  thrust  [Ref.  2].  No  damping  has  been  implemented.  The  thrust 
and  aerodynamic  forces  can  be  found  fairly  accurately  with  standard  test  and 
design  procedures  such  as  static  firings  to  obtain  thrust  versus  time  profiles  for 
the  engine  and  wind  tunnel  measurements  to  determine  the  aerodynamic  forces. 
Gravitational  forces  can  be  calculated  from  the  knowledge  of  the  missile's  posi- 
tion relative  to  the  earth.  The  mass  including  fuel  can  be  estimated  from  know- 
ledge of  the  missile's  weight  before  and  after  burnout  (from  the  measurements), 
and  by  using  a  mathematical  relationship  (often  linear)  for  the  decrease  in  missile 
mass  over  the  engine  burntime. 

In  this  work,  the  missile's  weight  is  assumed  constant  and  the  aerodynamic 
forces  are  zero.  The  deformations  resulting  from  the  structural  fiexibihty  have 
been  assumed  small  and  small  control  angle  assumptions  are  used  in  the  rigid- 
body  controller  design. 

The  primary  purpose  of  this  study  is  to  complete  the  simulation  of  the  flexible 
missile  in  2-D  motion,  where  the  bending  effect  is  considered  to  be  the  only  flex- 
ibility source  and  a  rigid-body  controller  as  developed  in  the  last  section  is  in- 
cluded. 

B.  SOLUTION  TECHNIQUE 

Many  numerical  integration  techniques  can  be  applied  in  the  solution  of  the 
equations  of  the  motion.  The  main  consideration  in  the  selection  of  the  inte- 
gration technique  is  the  size  of  the  time  step  necessary  to  integrate  the  equations 
of  motions  with  numerical  accuracy  and  stability. 

The  type  of  the  equations  of  motion  (2-20,  2-21)  in  this  research  permits  the 
application  of  the  Sequential  Integration  Method  [Ref.  6].  The  linear  equations 
of  small  motions  are  integrated  implicitly  and  the  large  motion  solutions  are  ob- 
tained using  explicit  integration  method.    The  implicit  methods  are  effective  for 
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linear  systems  with  high  frequencies  and  the  expHcit  methods  are  effective  for 
solving  nonlinear  systems  with  low  frequencies.  The  implicit  method  is  especially 
effective  for  hnear  systems  having  a  wide  range  of  frequencies  of  which  only  the 
lower  frequencies  are  excited.  The  size  of  the  time  step  need  only  be  chosen  to 
make  the  solution  of  the  excited  modes  sufficiently' accurate.  Explicit  methods 
are  effective  for  large  scale  systems  with  low  frequencies.  Because  of  the  low  fre- 
quencies, the  size  of  the  time  step  that  we  choose  for  stabihty  need  not  be  small. 
In  addition,  the  explicit  method  does  not  need'  iterative  procedures,  which  are 
time  consuming  for  non-linear  systems. 

C.  THE  COMPUTER  SIMULATION  CODE 

A  high  level  computer  language,  i.e.,  MATLAB,  was  chosen  to  simulate  the 
missile  system.  The  MATLAB  was  designed  for  matrix  operations.  The  simu- 
lation code  was  developed  with  modular  MATLAB  routines. 

The  simulation  code  can  be  divided  into  three  levels  :  LEVEL  1  (an  overview) 
separates  the  code  into  primary  portions  of  INITIALIZATION,  PLANT  DE- 
SCRIPTION, INTEGRATION,  SYSTEM  CONTROL  and  OUTPUT.  LEVEL 
2  facilitates  several  subroutines  for  LEVEL  1.  A  listing  of  MATLAB  routines 
required  for  manipulation  in  LEVEL  2  is  placed  in  LEVEL  3. 

D.  SIMULATION  RESULTS 

The  computer  simulation  was  performed  with  variable  parameters  which  de- 
termine missile  speed  and  controller  bandwidth.  These  parameters  include  force 
T,  K^  and  Kp. 

The  simulation  outputs  will  be  presented  in  large  motions,  small  motions  and 
control  angle.  The  large  motions  are  Xq,  Yq  and  6  and  the  small  motions  are  v(0) 
and  0(0).  The  initial  condition  for  all  runs  was  the  pitch  angle  of  45*.  The  sim- 
ulation work  was  divided  into  two  areas.  First  a  simulation  was  performed  for 
the  rigid  missile  using  the  rigid-body  controller.  These  results  were  used  as  a 
baseline  for  comparison  with  the  results  of  the  flexible  missile  using  the  rigid-body 
controller  system.  Second  a  simulation  was  performed  for  the  flexible  missile  and 
rigid-body  controller  system.  Only  the  trajectory  control  will  be  discussed  in  the 
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analyze  of  the  dynamics  of  the  flexible  missile.   The  desired  trajectory  was  speci- 
fied as 


e  = 


45°  yvhen  0.1  <r<0.l 

45"-  112.5r    when  0.1  ^r<0.4 

0  when  r  ^  0.4 


The  control  angle  is  assumed  to  be  limited  to  +/-10  degrees  (small  angle). 
This  restriction  puts  a  saturation  line  to  the  control  input. 

Fig. 5  represents  the  graphical  results  of  the  desired  and  actual  trajectory  for 
the  rigid  missile  and  rigid-body  controller.  The  force  (T)  and  controller  band- 
width (ct)„)  were  30000  N  and  3  rad/'s  respectively.  The  dark  black  line  and 
dashed  line  represent  the  desired  and  actual  trajectory,  respectively.  The  control 
angle  {b)  is  shown  in  Fig. 6  and  Figures  7-8  present  the  base  positions  {X^Y^  in 
large  motion. 

After  presenting  results  for  the  rigid  missile  with  rigid-body  controller,  the 
dynamic  behavior  of  the  flexible  missile  with  rigid-body  controller  will  be  studied 
next.  The  force  (T)  and  controller  bandwidth  (ct>„)  again  were  30000  N  and  3 
rad/s,  respectively.  Fig.9  shows  the  desired  and  actual  trajectory.  The  control 
angle  is  presented  in  Fig.  10.  After  O.l  sec,  the  missile  initially  needs  a  control 
angle  which  is  less  than  10  degrees.  The  effects  of  small  motion  can  be  seen  on 
control  angle  clearly.  The  base  deflection  and  slope  are  shown  in  Figures  11-12. 
There  is  no  small  motion  effects  on  flexible  missile  between  0-0.1  sec  because  of 
zero  control  angle  and  no  force  components.  After  0.1  sec,  the  small  motions  are 
excited  and  have  an  amplitude  of  10"'.  Figures  13-14  present  the  base  positions 
of  the  flexible  missile  in  large  motion.  The  large  motion  behaves  like  the  rigid- 
body  motion,  which  implies  that  the  coupling  effect  between  the  large  motion  and 
small  motion  is  small.  The  dynamic  behavior  of  the  flexible  missile  is  dominant 
by  the  large  motion  in  this  case  since  the  bandwidth  of  the  controller  is  low  and 
small  motion  is  not  significantly  excited. 
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Figure  5.      The  desired  and  actual  trajectory  for  the  rigid  missile  nitli  rigid-body 
controller  (T  =  30000  N,  a)„  =  3  rad/s) 
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Figure  6.      The  control  angle  for  the  rigid  missile  »ith  rigid-body  controller  (T 
30000  N,  aj„  =  3  rad/s) 
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Figure  7.      The  large  motion  position  X  for  the  rigid  missile  with  rigid-body  controller 

(T  =  30000  N,  co„  =  3  rad/s) 
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Figure  8.      The  large  motion  position  Y  for  the  rigid  missile  nith  rigid-body  controller 
(T  =  30000  N,  co„  =  3  rad/s) 
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Figure  9.      The  desired  and  actual  trajectory  for  the  flexible  missile  with  rigid-body 
controller  (T  =  30000  N,  a>„  =  3  rad/s) 
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rigure  10.      The  control  angle  for  the  flexible  missile  nith  rigid-body  controller  (T 

=  30000  N,  w„  =  3  rad/s) 
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Figure  1 1.      The  small  motion  position  v(0)  for  the  flexible  missile  with  rigid-body 
controller  (T  =  30000  N,w„  =  3  rad/s) 
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Figure  12.      The  small  motion  position  plii(O)  for  tlie  flexible  missile  >vith  rigid-body 
controller  (T  =  30000  N,  w„  =  3  rad/s) 
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Figure  13.      The  large  motion  position  X  for  the  flexible  missile  with  rigid-body  con- 
troller (T  =  30000  N,  a;„  =  3  rad/s) 
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Figure  14.      The  large  motion  position  Y  for  tlie  flexible  missile  nith  rigid-body  con- 
troller (T  =  30000  N,  co„  =  3  rad/s) 
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Two  kinds  of  tests  were  performed  to  further  the  study  of  the  dynamic  be- 
havior of  the  flexible  missile. 

The  first  test  was  to  study  dynamics  due  to  the  increase  of  the  control  system 
bandwidth  (a;„).  The  desired  and  actual  trajectory  of  the  fiexible  missile  with 
rigid-body  controller  without  saturation  on  the  control  angle  is  presented  in 
Fig.  15.  The  force  (T)  is  still  30000  N  while  the  controller  bandwidth  (co„)  is  in- 
creased to  300  rad/'s.  As  the  controller  bandwidth  is  increased,  the  gap  between 
controller  bandwidth  and  system  natural  frequency  is  smaller.  The  system  will 
be  unstable  at  controller  bandwidths  close  to  the  system  fundamental  frequency 
which  can  be  as  high  as  270  rad/s  in  this  case.  Note  that  the  fundamental  fre- 
quency of  the  simple  beam  with  free-free  boundary  conditions  was  calculated  as 
cl)„,  =209.2053  rad/s,  and  the  missile's  fundamental  frequency  will  be  increased 
due  to  the  coupling  between  large  and  small  motion.  The  control  dynamics  in- 
terfere with  the  structural  dynamics.  The  higher  control  frequency  with  the 
bandwidth  of  300  rad/s  significantly  excites  the  small  motion  of  the  structure. 
The  graphs  of  control  angle  and  the  small  motions  and  large  motions  of  the  mis- 
sile base  (Figures  16-17-18-19-20)  show  the  unstable  state. 

The  bandwidth  of  the  controller  must  be  set  much  lower  than  the  funda- 
mental frequency  of  the  missile  in  order  to  use  the  rigid-body  control.  This  implies 
that  the  pitch-angle  response  speed,  which  is  determined  by  the  control  band- 
width is  limited.  To  keep  the  response  speed,  the  missile  structure  must  be  de- 
signed sufficiently  rigid  to  possess  a  high  fundamental  frequency.  On  the  other 
hand,  the  missile  payload  will  affect  the  natural  frequency  of  the  missile  structure 
with  the  heavier  the  payload,  the  lower  the  fundamental  frequency  of  the  missile. 
Therefore,  the  payload  must  be  limited  to  achieve  high-speed  control  response. 
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Figure  15.      The  desired  and  actual  trajectory  for  the  flexible  missile  uith  rigid-body 
controller  without  saturation  on  the  control  angle  (T   =   30000  N,  (0„ 
.  =  300  rad/s) 
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Figure  16.  Tlfe  control  angle  for  the  flexible  missile  with  rigid-body  controller 
without  saturation  on  the  control  angle  (T  =  30000  N,  cu,  =  300 
rad/s) 
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Figure  17.  The  small  motion  position  v(0)  for  tlie  flexible  missile  rigid-body  con- 
troller uithout  saturation  on  the  control  angle  (T  =  30000  N,  (x)„  = 
300  rad/s) 
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Figure  18.  The  small  motion  position  phi(O)  for  the  flexible  missile  rigid-body  con- 
troller without  saturation  on  the  control  angle  (T  =  30000  N,  cu„  = 
300  rad/s) 
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Figure  19.  The  large  motion  position  X  for  the  flexible  missile  rigid-body  controller 
without  saturation  on  the  control  angle  (T  =  30000  N,  a)„  =  300 
rad/s) 
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Figure  20.  The  large  motion  position  Y  for  the  flexible  missile  rigid-body  controller 
without  saturation  on  the  control  angle  (T  =  30000  N,  co„  =  300 
rad/s) 
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Fig. 21  presents  the  desired  and  actual  trajectory  of  the  flexible  missile  using 
rigid-body  controller  with  a  saturation  ( +  /-10  Degrees)  on  the  control  angle.  The 
force  (T)  and  controller  bandwidth  (coj  are  30000  N  and  200  rad/s,  respectively. 
A  better  trajectory  tracking  is  obtained,  however  the  actual  trajectory  is  oscillat- 
ing about  the  desired  trajectory  which  is  mainly  due  to  the  switching  of  the  con- 
trol angle  between  the  saturation  lines.  Fig. 22  shows  the  control  angle,  and 
Figures  23-24  show  the  small  motion  displacement  and  slope  of  the  base  of  the 
flexible  missile.  It  was  observed  that  the  switching  of  the  control  angle  between 
+  /-10  degrees  causes  increased  excitation  of  small  motions.  Figures  25-26  pre- 
sents the  position  of  missile  base  in  large  motion  where  the  large  motion  is  no 
longer  dominated  by  the  rigid-body  motion. 

The  second  test  was  to  study  the  dynamics  of  the  flexible  missile  due  to  an 
increase  of  the  missile  speed.  To  increase  the  speed  of  the  missile,  the  thrust  force 
was  increased.  Fig. 27  presents  the  actual  and  desired  trajectory  of  flexible  mis- 
sile. The  force  was  60000  N  and  controller  bandwidth  was  3  rad/s.  As  seen  from 
Fig. 27,  the  pitch-angle  response  is  determined  directly  from  the  bandwidth  of  the 
controller.  The  increased  missile  speed  does  not  change  the  pattern  of  the  pitch- 
angle  response.  The  control  angle  of  the  flexible  missile  is  shown  in  Fig. 28.  The 
control  angle  is  affected  by  the  increased  speed.  When  the  speed  is  increased,  the 
control  angle  gets  smaller  because  a  smaller  control  angle  generates  a  sufficient 
correction  to  the  pitch  angle.  The  increased  missile  speed  has  little  effect  on  small 
motions  (Figures  29-30).  Figures  31-32  show  the  large  motions. 
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Figure  21.  The  desired  and  actual  trajectory  for  the  flexible  missile  using  rigid-body 
controller  with  saturation  on  the  control  angle  (T  =  30000  N,  co„  = 
200rad/s) 
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Figure  22.      The  control  angle  for  the  flexible  missile  using  rigid-body  controller  nilh 
saturation  on  the  control  angle  (T  =^  30000  N,  a;„  =  200  rad/'s) 
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Figure  23.  The  small  motion  position  v(0)  for  the  flexible  missile  using  rigid-body 
controller  >\ith  saturation  on  the  control  angle  (T  =  30000  N,  co„  = 
200  rad/s) 
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Figure  24.  The  small  motion  position  phi(O)  for  the  flexible  missile  using  rigid-body 
controller  with  saturation  on  the  control  angle  (T  =  30000  N,  a>„  = 
200  rad/s) 
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Figure  25.  The  large  motion  position  X  for  the  tlexibie  missile  using  rigid-body 
controller  nitli  saturation  on  the  control  angle  (T  =  30000  N,  0}„  = 
200  rad/s) 
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Figure  26.  The  large  motion  position  Y  for  the  flexible  missile  using  rigid-body 
controller  with  saturation  on  the  control  angle  (T  =  30000  N,  a)„ 
=  200  rad/s) 
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Figure  27.      The  desired  and  actual  trajectory  for  the  flexible  missile  using  rigid-body 
controller  with  increased  speed  (T  =  60000  N,  co„  =  3  rad/s) 
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Figure  28.      The  control  angle  for  the  flexible  missile  using  rigid-body  controller  uith 
increased  speed  (T  =  60000  N,  a>„  =  3  rad/s) 
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rigure  29.      The  small  motion  position  v(0)  for  the  flexible  missile  using  rigid-body 
controller  with  increased  speed.   (T  =  60000  N,  cu„  =  3  rad/s) 
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Figure  30.      The  small  motion  position  plii(O)  for  the  flexible  missile  using  rigid-body 
controller  nitli  Increased  speed  (T  =  60000  N,  a)„  =  3  rad/s) 
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Figure  31. 


The  large  motion  position  X  for  the  flexible  missile  using  rigid-body 
controller  with  Increased  speed  (T  =  60000  N,  aj„  =  3  rad/s) 
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Figure  32.       The  large  motion  position  Y  for  the  flexible  missile  using  rigid-body 
controller  with  increased  speed  (T  =  60000  N,  a)„  =  3  rad/s) 
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V.     SUMMARY 

A.  CONCLUSIONS 

The  development  and  computer  simulation  is  presented  for  a  bending  flexible 
missile  with  a  rigid-body  controller  system  moving  in  a  2-D  coordinate  frame.  In 
this  research,  the  dynamic  model  has  been  developed  through  the  Equivalent 
Rigid  Link  System  utilizing  Lagrangian  dynamics  to  obtain  a  type  of  system 
equations  of  motion  suited  for  Sequential  Integration  Method  that  integrates 
large  motion  explicitly  and  small  motion  implicitly.  The  spatial  finite  element 
discretization  of  missile  structure  and  the  application  of  truncated  natural  modal 
responses  provide  an  approximate  solution. 

The  analysis  and  simulation  were  performed  to  understand  the  dynamic  be- 
havior  of  a  Hexible  missile  using  a  rigid-body  controller.  It  was  found  that  the 
controller  bandwidth  must  be  much  lower  than  the  fundamental  frequency  of  the 
missile  in  order  to  use  the  rigid-body  controller.  The  payload  will  affect  the  na- 
tural frequency  of  the  missile  structure  i.e,  when  the  payload  is  increased,  the 
system  fundamental  frequency  will  be  decreased.  The  payload  must  then  be  lim- 
ited to  achieve  high-speed  response.  In  order  to  increase  the  payload  and  main- 
tain high-speed  control  response,  a  flexible-body  controller  is  needed. 

B.  RECOMMENDATIONS 

Areas  which  remain  to  be  investigated  include  : 

l.Add  the  payload  and  aerodynamic  effects  to  the  model. 

2.Design   and   study   the  dynamic   behavior  of  a   flexible   missile  with   a 

flexible-body  controller  in  2-D  motion. 
3. Build  a  Hexible  missile  in  a  laboratory  scale  and  obtain  experimental  data. 
4. Design  and  simulate  a  control  system  for  flexible  missiles  in  3-D  motion. 
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APPENDIX  A.     DERIVATION  OF  THE  MODE  SHAPE  RESPONSE 
MATRIX  FOR  THE  FLEXIBLE  MISSILE 

In  this  study,  the  natural  mode  shape  functions  of  a  beam  are  used  to  repre- 
sent the  flexural  motion  of  the  flexible  missile.  Only  the  first  two  mode  shapes  are 
used.  The  flexible  missile  is  modeled  as  a  continuous  Euler-Bernoulli  free-free 
beam,  neglecting  shear  deformation  and  rotary  inertia  effects. 

The  bar  (Fig. 33)  has  system  parameters: 


Figure  33.      The  Bar  in  Flexure 

y(x,t)  =  Transverse  displacement  at  any  point  x  and  time  t 

f(x,t)  =  Transverse  force  per  unit  length 

m(x)  =  The  mass  per  unit  length 

EI(x)  =  Flexural  rigidity 

E       =  Young's  modulus  of  elasticity 

I(x)    =The  cross-sectional  area  moment  of  inertia 

Q(x,t)  =  Shearing  force 

M(x,t)  =  Bending  moment 
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Fig.34  shows  the  free  body  diagram  corresponding  to  a  bar  element  of  length 


dx. 


Figure  34.      Free  Body  Diagram  Corresponding  to  a  Bar  Element 


The  force  balance  of  the  free  body  is 


dQixJ)  d  y(x,t) 

iQ{x,t)  +  -^-^  ^]  -  Q{x,t)  +Ax.t)dx  =  m{x)dx     ^\^ 
dx  ^t^ 


{A-\) 


The  moment  equation  of  motion  about  the  axis  normal  to  x  and  y,  ignoring 
the  inertia  torque  associated  with  rotation, 
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dMix,t)  dQ(xj) 

lM{x,t)  +  — j^  dx^  -  M{x,t)  +  lQ{x,t)  +     ^j^   '  dxyx 

+  J{x,t)dx-^  =  0  {A -2) 


Canceling  appropriate  terms,  ignoring  terms  involving  second  order  terms  in 
dx  and  combining  Eqs.  (A-1)  and  (A-2), 

d\M{x,t)    ,   ^     ^      ,^^     ^  d^y{x,t) 


^  +A^,t)  =  M{x,t)     ■"';'  [A  -  3) 

ex  ct 

Eq.  (A-3)  relates  the  bending  moment  M(x,t),  transverse  force  f(x,t)  and 
bending  displacement  y(x,t). 

The  relation  between  the  bending  moment  and  the  bending  deformation  is 

M{x,t)  =  EI{x,t)  —-4^  [A  -  4) 

dx 

Inserting  Eq.  (A-3)  into  Eq.  (A-4),  we  obtain  the  differential  equation  for  the 
flexural  vibration  of  a  bar, 

r-  (£/(x) r— )  +yi-x,r)  =  M{x,t) —  {A  -  5) 

hx^  dx^  dr 

The  bending  moment  and  shearing  force  of  the  free  ends  (x  =  0,  x  =  L)  are 
zero. 

dx^  dx^ 

Eqs.  (A-6)  and  (A-7)  are  called  natural  boundary  conditions. 
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Considering  the  free  vibration  characterized  by  f(x,t)=0,  using  separation  of 
variables  method  and  simplifying  Eq.  (A-5), 

"^     \_EI{x)  m^  ]  =  w^mix)  Y{x)  {A  -  8) 


dx'^  dx^ 

Simplifying  the  Eq.  (A-8),  for  EI(x)=  constant 

dx 

where  /?*  =  ^  ^ 


£/(jc) 

The  boundary  conditions  require  that, 
At  x  =  0  (base) 


d'Y{x)  ^       _ 


2       1^=0 

ax 

^Uo^O  (.-.0, 

dx 

At  x  =  L  (tip) 


etc'      '"^ 


'^'^'^•.,.-0  (A-U) 


The  general  solution  of  Eq.  (A-9), 
Y{x)  =  Ci  sin(/?x)  +  C2  cos(/?;c)  +  C3  sinh(/?x)  +  C4  cosh(/?x)  {A  -  12) 

Taking  the  derivatives  of  Eq.  (A- 12),  and  substuting  boundary  conditions, 
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Y{x)  =  Ci(  sin  px  +  sinh  {]x)  +  C2(  cos  fJx  +  cosh  px)  (^  -  13) 

Solving  Eq.  (A- 13)  yields, 
cos /?L  cosh /IL  =  1  {A-\4) 

The  first  five  consecufive  roots  of  this  equation  are; 

p,  =  0.0 

/?2  =  4.712388 

/?3  =  7.853981 

p,=  10.995574 

Ps=  14.137166 

Eq.  (A- 13)  is  now  written  in  the  following  form, 

Y^{x)  =  C^(  cos  p^  +  cosh  p^)  +  ( sin  p^  +  sinh  p^)     r  =  \,2  (^-15) 

where: 

sin  p.L  —  sinh  p^L 
'~  -  cos  p^L  +  cosh  p,L  ^ 

The  transverse  displacement  v(x)  and  slope  <l)(x)  can  be  represented  in  the 
following  forms  respectively, 

v(.x)=  l^a.YXx)  {A-\l) 

cDW  =  -g-  (^-18) 

v(.x)  =  fli(Ci(  cos  P^x  +  cosh  P^^x)  +  ( sin  P^x  +  sinh  P\x)) 

+  «2(Q(  cos  /J2^  +  cosh  /?2.x)  +  ( sin  /?2^  +  sinh  PiA)  (^  —  19) 

cl)(x)  =  a^iC^PyX  +  sinh  j^^x)  +  P^  cos  /?ix  +  cosh  P^x)) 

+  ^2(<^2i^2(  -  sin  /?2-x:  +  sinh  /^jjc)  +  ( cos  P-^x  +  cosh  /Jjx))  (/I  -  20) 
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Substituting  the  boundary  conditions  into  shape  function  gives, 

v(0)  =  la^C^  +  2^2^  (^  -  21) 

0(0)  =  2^1/?!  +  2a2^2  (^  -  22) 

The  modal  amplitudes  a,  and  a2', 

^2  =  -  -^  v(0)  +  -^  cD(O)  {A  -  24) 

V  =  tzv(O)  +  Z>(D(0)  .__  (^-25) 

a  =  FJ^{x)-\-F,Y2{x)  {A -16) 

b  =  F2Y,{x)  +  F^Y2{x)  [A -21) 


-2 


2<-i        crE 


-c 


2 


F2  =  -r^  {A-29) 

F,  =  -^  ■  (^-30) 


1 

E 


F^  =  -^  {A-3\) 


E  =  2p2-2^  {A-32) 


Substituting   Eqs.  (A-15),  (A-16),  (A-26),  (A-27),  (A-28),  (A-29),  (A-30), 
(A-31),  (A-32)  into  Eq.  (A-25)  yields  an  expression  for  the  transverse  displace- 
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ment  of  a  flexible  missile  as  a  function  of  the  missile  base  nodal  displacements, 
v(0)  and  0(0), 

v{x,t)  =  [F^{Ci{  cos  /?ix  +  cosh  p^x)  +  ( sin  PiX  +  sinh  Pix)) 
+  F3(C2(  cos  /?2^  +  cosh  /?2-x)  +  ( sin  /?2-X"  +  sinh  /?2-^))]^(0) 
+  [FiiCii  cos  /?ijc  +  cosh  Pix)  +  ( sin  P^x  +  sinh  j^i^)) 


+  F4(C2(  cos  /?2X  +  cosh  /?2^)  +  ( sin  /?2^  +  sinh  /?2^))]^(0) 


(/4  -  33) 


This  expression  is  differentiated  twice  to  obtain  v"{x),  which  is  necessary  for 
the  calculation  of  the  potential  energy  due  to  deformation  and  theoretical  strain, 

v"(;c,0  =  {FiiCiPli  -  cos  P^x  +  cosh  P^)  +  pl{  -  sin  p^x  +  sinh  P^x)) 
+  ^3(C2/^2(  -  COS  P2X  +  cosh  P2)  +  Pli  -  sin  P2X  +  sinh  /J2-v)))v(0) 
-I-  [FiiCiPli  -  cos  Pix  +  cosh  Pi)  +  pI{  -  sin  P^x  +  sinh  i^^x)) 


+  ^4(Q/^2(  -  cos  /?2X  +  cosh  P2)  +  /?2(  "  sln  P2X  +  sinh  /?2X)))O(0) 


{A  -  35) 


Now  substitution  of  v(x)  into  the  3x1  deformation  vector,  d  ,  yields  the  3x2 
shape  function  matrix,  </),  and  the  2x1  nodal  displacement  vector,  U, 


d  =  (pv- 


'  0 

"0  o' 

'  v(0)  ■ 

0 

— 

0  0 

_CD(0)_ 

_v{x)_ 

a  b 

{A  -  36) 


The  shape  function  matrix  is  now  in  a  form  convenient  for  computer  coding. 
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APPENDIX  B.     DERIVATION  OF  THE  EQUATIONS  OF  MOTION  FOR 

THE  FLEXIBLE  MISSILE 

The  ERLS  is  used  to  separate  the  flexible  missile's  motion  into  a  small  motion 
and  a  large  motion.  The  large  motion  generalized  coordinates  are  defined  as  X^, 
Yq,  theta  (O)and  the  small  motion  generalized  coordinates  as  U.  The  Lagrangian 
Dynamics  are  very  helpful  in  the  derivation  of  the  equations  of  motion  of  com- 
plex systems.  U  represents  a  column  vector  which  has  two  elements.  First  one  is 
v(0)  which  is  the  nodal  displacement  of  the  flexible  missile's  base  and  ^{Q)  which 
is  the  slope  of  the  flexible  missile's  base.  We  can  represent  the  vector  of  the  gen- 
eralized coordinates  by  Eq.  (B-1), 

1  =  1    X^    Y^    e    y  t/  =  [v(0)    (D(0)    y  (5-1) 

These  generalized  coordinates  will  be  used  in  the  application  of  the  Lagrange 
equations  to  the  equations  of  motion  for  both  small  and  large  motion  coordinate 
stems.   The  Lagrangian  equations  are  defined  in  Eqs.  (B-2)  and  (B-3). 

d  r-  cKE  -1  cKE    ,    oPE       r-   /•   i  -^  -}\  (r>     '^\ 

-17\.—-—^-~J7~^~Tr^^^i    0=^'23)  (5-2) 


d   r  dKE  -,  dKE    ,  dPE       ^  ro   l^ 

[ —  J IT-  +  — =^  =  fu  (5  -  3) 


"^^         -.  dU         dU 

dU 

where  i^,  is  a  component  of  ^ 

Using  ERLS  we  can  define  the  global  position  of  the  base  position  in  terms 
of  the  local  position  vector  (r),  a  deformation  vector  (d),  and  a  coordinate  trans- 
formation matrix  (W)  with  Eq.  (B-4), 

R  =  W{T  +  d)  {B-  4) 
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Continuing  in  the  development  of  the  equations  of  motions,  we  determine  the 
derivative  with  respect  to  time  for  global  positions  in  order  to  obtain  the  kinetic 
energy  expressions. 

R  =  W{T  +  d)+  Wd  {B  -  5) 

The  kinetic  and  potential  energy  expresions  follow.  In  this  development  we 
have  seperated  the  large  and  small  motion  energy  contributions  and  will  present 
them  seperately. 

KE  =  Yi^^^dm  {B-6) 

PE  =  Y  jt/^r^CrC^i/jc  -  \R^gdm  [B  -  7) 

Redefining  d  in  terms  of  the  shape  matrix  (f)  (derivation  of  d  is  presented  in 
Appendix  A  ),  we  can  rewrite  Eq.  (B-6)  as 

KE  =  y  j{[lV{T  +  d)+  WdflW{7  +  ?)  +  lVd^)dm  {B  -  8) 


+  U^(1)^W'^W4>U  +  2U'^4)'^W'^W4)U  +  u'^4)'^W^W4)U^)dm  [B  -  9) 

Continuing  in  the  development  of  the  equations  of  motion,  we  first  express 
the  derivative  of  the  kinetic  energy  with  respect  to  the  time  rate  of  change  of  the 
large  motion  coordinate  ^,  Eq.  (B-13),  and  then  determine  the  time  rate  of  change 
of  this  expression,  Eq.  (B-14).   Since 
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^  = 


■^^o' 

"c^r 

^'o 

= 

^2 

_  e_ 

^3 

(5-10) 


VV,= 


dW       dW 


n 


,2 


IV.  = 


dX 

d^, 

dlV 

div 

dY 

d^2 

dlV 

dlV 

dd 


d^, 


=  Partial  derivative  w.r.t.  X 


=  Partial  derivative  w.r.t.  Y 


=  Partial  derivative  w.r.t.  d 


IF,  =  Partial  derivative  of  W  w.r.t.  X,  Y,  0   i  =1,2,3 


— : —  =  -r-  \{r    : —  Wr  +  r   IV   — : —  r  +  2r    : —  WoU 


d^i 


5^. 


,-T,i,T  dlV    ,-  ,   ^-T  dlV^ 


W({)U  +  U'cf) 


T_,T  d\v'   ,;.^r: 


5q 


W4)U 


■jrT 


W(f)U)dm 


(5-11) 


Putting  Eq.  (B-12)  into  Eq.  (B-11)  and  simplifying 


dW       dlV 


d^i        ^^i 


=  W, 


dii  ^^i 


(5-12) 


dKE      r^-TTwrT, 


TtifT, 


^^i 


^  =  J(  PwJwT  +  rwJW(pu  +  7^w^Wi<i)U 


+  T'^lvJW(j)U  +  u'^(i)'^wJW(t)U  +  U'^4>'^W]W(1)U  )  dm 


(5-13) 
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+  T'^wJw4)U  +  T'^wJiV(f)U  +  T'^w^Wj(f)U 


+  7'^w'^w^^<i)U  4-  T'^w'^Wj<pu  +  T'^iVjiV(j)U 


+  r^wJW(t>U  +  7'^WJW4)U  +  u'^(l)^wJlV(l)U 


+  U'^(1)^WJ^IV(J)U  +  U^(l)'^lvJW(pU  +  U^(f)'^w]W(j)U 


+  u'^(i)^wJw4)U  +  u^(f)^w^iW(i)U  +  u^(i)^ivJiv4)U 


+  U'^(J)'^IVJW4)U  )dm  [B-  14) 

Finally  we  express  the  derivative  of  the  kinetic  energy  with  respect  to  the 
large  motion  position  as  c,. 

dKE  {/-TjirTjir^      ,     --TrirTrir,rr     ,     -^T  jirT 


d^i 


=  j(  T^w^wT  +  T^w;^W(j)U  +  T^w^w^^^u 


+  7'^wlw4)U  +  T^vv'^w^(pu  +  u'^(i)'^w'^^W(f)U 


+  u'^(})'^wlW(f)U  +  U'^(1)'^W^W,(1)U  +  U^4>'^W'^W^^4>U  )dm  [B  -  \5) 
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-J-  ( -^^ )  -  -^jf-  =  j(  7^w]w7  +  7^w]w4)U  +  2T^ivJw<i)U 


dti  ^C/ 


+  r'^W^Wj<t)U  +  7'^wJlV4)U  +  u'^(t>'^wJW(f)U 


+  2U'^(l)'^wJlV4)U  +  u'^(j)'^wJW(t>U  )dm  (5-16) 

Completing  our  development,  the  expressions  for  the  potential  energy  are 
presented  organized  similarly  to  the  previous  material, 

i^  =  -  j(f  +  (pufwjgdm  =  -  \7^wjgdm  -  \u'^(l)'^ivjgdm  (5-17) 

The  small  motion  equations  are  handled  in  a  similar  way,  and  are  easy  to 
derive  than  large  motion. 


dKE        [/ iTjj.Tjir-   ,    xTti/Tt;^,  w   ,     ,T-,j,T 


=  \{(^^W^W7  +  (l)W^W4>U  +  4)W^W(^U)dm  (5-18) 


dV 


tT 


^Ml.  =  \[(l)'^w'rw7  +  4)'^W'^W(I)U  +  <l)'^W^W<t)U)dm  {B  -  19) 


-J-  ( -^4^ )  =  \{(^'^w'^w7  +  (^'^w'^w7  +  <f)'^iv'^W(})U  +  (i>'^w'^w4>u 


dU' 


+2(f)'^w'^W(t)U  +  (l)'^lv'^W(t)U  +  <i)'^w'^W4>U  )dm  {B  -  20) 


du'^ 


dW 
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+  2(l)'^W^W(l)U  +  <j)'^\v'^W4>U  )dm  (5-21) 


Potential  energy, 
^f-  =  Iv'^CYUdx  -  jcp'^W^gdm  {B  -  22) 

We  can  simplify  these  expressions  furtlier  by  the  substation  of  the  second 
time  derivative  of  the  coordinate  transformation  matrix  W,  The  contributions  to 
the  coriolis  and  centripetal  forces  represented  by  IV,  ,  termed  the  residual  accel- 
erations. The  contribution  to  the  general  forces  is  represented  by  IVj . 

W=W,+  Y^Wjti  (5-23) 

For  large  motion, 


y=l  '  '  '    J=l 


+  2U'^(})'^lvJlV4)U  +  Ur4>^wJw4)U  -  y^wjg  -g^W^^(l)U)dm  =  F^         {B  -  24) 


For  small  motion, 


J=i  j=i     '  '        , 
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+  <p^w'^w<i)U  +  r^cru  -  (jy'^w^g  =  f^ 


[B  -  25) 


Collecting  terms  and  arranging  the  coefficients  the  two  Lagrange  equations 
can  be  written  as  the  equations  of  motion  for  the  large  and  small  generalized  co- 
ordinates. It  is  these  equations  which  must  be  solved  by  computer  simulation 
code. 

For  the  large  motion, 


3 

z 


\u'^(f)'^lvJlVj(t>Udm  + 
\U^(j)'^lvJlV/dm 


'^J  + 


\r'^lvJW(j)dm  + 
\u'^(f)'^lvJW(j)dm 


U  = 


F,j  -  \r'^wJw;Pdm  -  IT'^wJw^cpUdm  -  2\7'^lvJlV(t)Udm 
-  \U^(i)'^wJw;Pdm  -  lu'^<p^lvJw^(f)Udm  -2\U^(t)'^wJW(f)Udm 

+  \7'^ivjgdm  +  lg'^lV^^(j)Udm 


[B  -  26) 


For  the  small  motion, 


Z  \.\<l)^W^W/dm  +  l<p'^W^W^4)Udm^ij  +  l\(p'^w'^W(t)dm^U  + 


lJ2<j>^W^W(f)dm']U  +  [l(t)^lV,(t)dm  +  \rcrdx}U 


TrrrTjirr 


\(})^W^gdm  -  \4)^W^W/dm  +  F„ 


[B  -  27) 
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By  defining 


!^TrrrT 


TiTjirT, 


^qq{iJ)  =  \r^JlV/dm  +  i7^wJWj(j)Udm  +  lU^(J)^wJWj<pUdm 


Mqn  = 


^Kq  =  ^Un 


+  \W<i)^w]w/dm       (/=  1,2,3)         (/■=  1,2,3) 


\7^wlW(^dm  +  \u'^(^'^wlW4)dm 
\T'^W^W(t)dm  +  \u'^(t>'^[vlW({)dm 
\'i^lV^^W(i)dm  +  \U^(j)'^wlW(j)dm 


T 


M^^  =  \4>^lV^W(j)dm 


G^  =  \2(j)^W^W4)dm 


T, 


K^  =  \(})^iv^(f)dm  +  \rcrdx 


[B  -  28) 


{B  -  29) 


{B  -  30) 
(5-31) 
{B  -  32) 
(5-33) 


^T,,,T,i; 


f-T,j,Trl 


F^.  =  F,j  -  \rlvJlv;Pdm  -  yw]w^(l)Udm  -  2\7^wJlV(})Udm 


'TiTjj.Tij 


T   iTrrrT, 


-  \U^4>^w]W/dm  -  \W(l)^w]w^4)Udm  -2\W(l)^w]\V4)Udm 


rrr 


+  W^vjgdm  +  \g^W^^Udm  [i  =  1 ,2,3) 


{B  -  34) 


F„  =  \4)^W^gdm  -  \4>^W^w;Pdm  +  F„ 


The  final  equations  of  motion  are  written  as 


{B  -  35) 


{B  -  36) 
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M,,  ^  +  M,„  U    =  F^ 


{B  -  37) 


M,,  ^  +  M,,  t/  +  G,  ^  +  /:,  t/    =  F„ 


(5  -  38) 
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APPENDIX  C.     DERIVATION  OF  THE  GENERALIZED  FORCES  OF 

THE  EQUATIONS  OF  THE  MOTION 

To  derive  the  generalized  forces  the  principle  of  virtual  work  is  used.  We  first 
write  the  transformation  between  global  coordinates  and  local  coordinates,  Eqs. 
(C-1)  and  (C-2)  respectively, 


1 

X 

Y 

r 

y 


1         0  0 

.Yo    cos(O)    sin(0) 
Yq    s\n{d)    cos(0) 


1 

X 

y 


1  0  0 

-^0  cos(O)  -  Fo  sin(0)     cos(0)     sin(O) 
.Yo  sin(0)  -  y-Q  cos(0)     -  sin(0)    cos(0) 


1 

X 

Y 


(C-1) 


(C-2) 


From  Eqs.  (C-1)  and  (C-2)  we  can  find  x  and  y. 
JC  =  -.Yo  cos(O)  -  Fo  sin(O)  +  A'cos(0)  +  F sin(0) 

y  =  X^  sin(a)  -  Y^  sin(0)  +  .Ysin(O)  +  F  cos(O) 


(C-3) 
(C-4) 


Taking  derivative  of  the  Eq.  (C-3)  and  (C-4)  doing  necessary  substations,  we 
can  find, 


hx  =  hXzQs{Q)-^bYsm{Q) 

^y=-<5A'sin(^)  +  ^rcos(0) 

Using  the  virtual  work  principle, 
bW^  =  {_  -Tcos{5  +  cD(0))v(0)]((50  +  50(0)) 


(C-5) 
(C-6) 


+  lTcos{8  +  O(0))]5x  -h  [rsin(5  +  <I)(0))]((5y  +  5v(0)) 


(C-7) 


Based  on  small  angle  assumptions, 
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cos(^  +  a>(0))  =  1 


(C  -  8) 


and 


sin(5  +  (l>(0))  =  <5  +  <I>(0) 


(C-9) 


Doing  the  necessary  operations  and  substitutions, 
6U)^  =  (Tcosid)-T5  sin(0)-rcD(0)  smie))5X  +  (rsin(0)  +  T5  cos(0)  + 

rcl)(0)  cos{d))5  r  +  (  -  Tv{0))5{d)  +  (7(5  -f  rcl)(0))5v(0) 


+  {-Tv{0)5<t>{0)) 


(C-9) 


From  Eq.  (C-10),  one  can  write  large  motion  generalized  forces  (Fj),  i.e.,  Eq. 
(C-I  I)  and  small  motion  generalized  forces  (FJ,  i.e.,  equation  (C-12), 


f, 


Tcos{e)  -  T5  sin(0)  -  70(0)  sin(0) 

7'sin(0)  +  TS  cos(6))  +7'a)(0)  cos(0) 

-rv(0) 


(C-11) 


/^u  = 


r^  +  7'<i>(0) 
-rv(o) 


(C-12) 
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APPENDIX  D.     COMPUTER  SIMULATION  CODE 


"4%%%%%%%%%%%%%%%%%%%%%%%%%%  DEFINITICNS  CF  VARIi\BLES  %%%%%%%%%%%%%%%%%%%%%%%% 

%  % 

%  aO-7  =The  coefficients  for  the  numerical  integration  % 

%  area  =Cross  sectional  area  of  the  missile  (m**2)  % 

%  betal  =Mxial   node  1  % 

%  beta2  =Modal  node  2  % 

%  count  =Constant  that  defines  the  eleirent  positions  in  vectors.  % 

%  dtheta  =Desired  pitch  angle (dgrs)  % 

%  dthetadot=Desired  angular  velocity (rd/s)  % 

dtheta2dDt=C)e3ired  angular  acceleration ( rd/s **2)  % 

delta  =Control  angle (dgrs)  % 

E  =Young's  nodules (pascal)  % 

%  epsilon  =Parameter  for  integration  coefficient  % 

%  fn  =Shv3ll  motion  force  coefficient  matrix  :   % 

%  flexible=Constant  that  does  the  flexible  part  on  and  off        '  % 

"6  fnn  =Befonned  small  motion  force  coefficient  matrix  % 

1  fq  =large  motion  force  vector  % 

%  ftime  =Finish  tiine  of  the  integration  % 

%  force  =Elxtemal  force  (Newton)  % 

%  gn  =Gyrosccplc  coefficient  matrix  % 

%  grav  =Gravitational  constant (m/s**2)  % 

%  grawec  =Gravity  matrix  ;    % 

%  h'  =Tlme  st^  used  in  the  integration  % 

%  Izz  =Mcinent  of  inertia  about  the  z-axis(m**4)         /  % 

%  iota  =Parameter  for  integration  coefficients  % 

l  kn  =Stifness  matrix  % 

t  k  =NuniDer  of  the  column  of  the  matrix  % 

%  Icp,  )cv  =The  position  and  velocity  feedback  gains  % 

%  1  =Munt5er  of  the  row  of  the  matrix  % 

%  1 force  =Large  motion  force  vector  % 

1;  IjtI  =The  length  of  the  missile  (meter)  % 

%  Imn  =Constant  (lm*n)  % 

%  mn  =Feformed  small  motion  inertia  and  coupling  coefficient  matrix   % 

%  mqn  =Coupling  inertia  coefficient  matrix  % 

%  mqq  =Large  motion  inertia  coefficient  matrix  % 

%  ma  =Constant  (mu*area)  % 

%  nphi  =Shape  function  matrix  % 

"h  rtphi2  =Second  time  derivative  the  shape  function  % 

%  mu  =^1as3  density  of  the  missile (kg/m**3)  % 

%  n  =NuiTt)er  of  the  Integration  % 

%  phi  =Small  motion  slope  position  % 

pliidot  =anall  motion  slope  velocity  % 

%  phl2dot  =ShBll  motion  slope  acceleration  % 

%  rigid  =Constant  that  does  the  controller  on  and  off  % 

"k  rlocat  =Local  position  vector  % 
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q. 


%  q              =Large  motion  generalized  position  vector  % 

%  qckit         =Large  motion  generalized  velocity  vector  % 

%  q2dot      =Large  motion  generalized  acceleration  vector  % 

%  stlf        =Stlfness  matrix  % 

%  sforce    =Sinall  motion  force  vector  % 

%  slope          =Trajectory  path  slcpe  % 

%  tenpflrst=Teiiporary  matrix  for  the  Sinpson's  rule  % 

%  tenplast=TeiTporary  matrix  for  the  Sinpson's  rule  % 

%  tenpl ...  % 

%  tenp6    =Tenporary  matrices  for  integration  % 

%  tenp7 ...  % 

%  tenplO   =Tenporary  matrices  for  large  and  small  motion  coefficients  % 

%  theta   =Large  motion  angular  position  (pitch  angle-dgrs)  % 

%  thetadot=Large  motion  ang\ilar  velocity  (rd/s)  % 

%  theta2dot=Large  motion  angular  acceleration (rd/s**2)  % 

%  time    =Sec  % 

%  u      =Small  motion  generalized  position  vector  % 

%  udot    =Small  motion  generalized  velocity  vector  % 

%  u2dot   =3Tiall  motion  generalized  acceleration  vector  % 

%  urn     =aTBll  motion  position  vector  % 

%  umdot   =Small  motion  velocity  vector  •.  % 

%  V      =Small  motion  deflection  position (m)  % 

%  \Kiot         =Sinall  motion  deflection  velocity (m/s)  % 

%  v2dot   =Sinall  motion  deflection  acceleration  (m/s**2)  % 

%  wresld  =Residual  acceleration  matrix  % 

%  wdot    =Time  derivative  of  the  transformation  matrix  % 

%  wksl    =The  derivative  of  the  transformation  matrix  with  respect  to  the  % 

%         large  motion  generalized  coordinates  % 

%  X      =Large  motion  x-dlrectlon  position (m)  % 

%  X      =Local  X  coordinate  % 

%  xcjot    =Large  motion  x-dlrection  velocity  (m/s)           ^/  % 

%  x2dot   =Large  motion  x-dlrection  acceleration (m/s**2)  % 

%  y      =Large  motion  y-dlrectlon  position (m)  % 

%  ydot    =Large  motion  y-directlon  velocity (m/s)  % 

%  y2dot   =Large  motion  y-dlrection  acceleration (m/s**2)  % 

%  mnnl,)cnl,  % 

%  )<n2,fnl,  % 

%  fn2,fn33,  % 

%  gnl,mc(ql,  % 

%  fql,fq2  =Tenporary  matrices  for  the  calculation  of  the  coefficients  % 

%  cl,c2,fl,  % 

%  f2,  f 3, f 4=Shape  function  constants  % 

%   m22,f22,  % 

%  h2,b2,A,  % 

%  B,F,terpl4,  % 

%  tenp20  =Tenporary  matrices  for  the  controller  % 

%  mll,ml2,  % 

%  m21,m22,  % 

%  fll,hl,bl=Constant3  for  the  controller  % 

%  tDlot    =Tline  plot  vector  % 
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%  xplot    "=Large  motion  x-dlrection  position  plot  vector  % 

%  yplot    =Large  motion  y-dlrection  position  plot  vector  % 

%  thetaplot=Large  motion  angular  position  plot  vector  % 

%  xdotplot  =Large  motion  x-direction  velocity  plot  vector  % 

%  ydotplot  =Large  motion  y-direction  velocity  plot  vector  % 

%  thetadotplot=Large  motion  angular  velocity  plot  vector  % 

%  vplot    =Shiall  motion  deflection  position  plot  vector  % 

%  phiplot  =anall  motion  slope  position  plot  vector  % 

%  vdotplot  =Shiall  motion  deflection  velocity  plot  vector  % 

%  phldotplot=Sniall  motion  slope  velocity  plot  vector  % 

%  Deltaplot  =Control  angle  plot  vector  % 

%  dthetaplot=Desired  trajectory  angular  position  plot  vector  % 

%  dthetadeg, deltadeg,  % 

%  thetadeg  =Angle3  converted  frcm  radian  to  degrees  % 

%  % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  MMN  PROGRAM  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  LEVEL  1  % 

%  This  is  the  controlling  program  for  flexible  missile.  It  defines  the  % 
%  variables,  determines  the  coefficients  for  explicit  integration  and  mode  % 
%  shape  function,  caj-ls  large  and  small  motion  coefficient  matrices  routines, % 
%  large  emd  small  motion  Integration  routines,  controller  routines,  output  % 
%  routines  % 

%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


(X, y,  theta,xdot,ydot,thetadot,  lm,betal,beta2,E,  Izz,  force,  . . . 
phi,  delta, V,  lota,epsllon, ftime, h,grav,ma, vdDt,phidot,n,  Imn, flexible, 
count, nphl,riphl2,stlf,gravvec,rlocal,x2cbt,y2dot,theta2dot,v2dot,  . . . 
phi2dot]  =  constantsl; 

(m22,  f 22,  h2, b2.  A,  F,  B,  dtheta2dot,  dthetadot, dtheta,  kv,  . . . 

kp,  rigid,  tenpl4,  slcpe,  terp20]  =constantsla; 
(q,qdot,q2dot,tG^p6,wresid,wdot,wksl,w,  . . . 

Iforce,  s force, um,umdot,fnn, gn,kn,mn, count  1, mil, ml2,m21,  fll,hl,  . . . 

bl  ]  =constants2a  (X,  y,  theta,  xdot,  ydot,  thetadot,  x2dot,  y2dot,  theta2dot)  ; 

(u,udot,u2dot,mnnl,knl,kn2, fnl, fn2, fn33,mnn, gnl,mqq,mqn, . . . 
fq,mqql,tenpl,tenp2,tenp3,tefip4,terp5,teitp7,mqnl,terp8,  fql,tenrp9,  . , . 
tetrplO,  fq2J  =  constants2b(v,Fiil,vdot,phldot,v2dot,phl2dot); 

( tplot, xplot, yplot, . . . 

thetaplot,xdotplot,ydotplot,th^tadotplot, vplot, phiplot,  vdotplot, . . . 

phidotplot, thetadeg,  deltaplot, dthetaplot,dthetadeg,deltadeg]  = 

constantsB ( f time,  h) ; 

[a0,al,a2,a3,a4,a5,a6,a7]  =  ccoef (epsllon,  lota,h)  ; 
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(cl,c2,  fl,  12,  f3,  f4]=Tiphicons  (ljn,betal,beta2)  ; 

for  time  =  O.Orhrftlme, 

tw,wckDt,wksl,wre3id,lforce,sforce,im,undot)   =  wmatderlwariant  ( . . . 
theta,  X,  y,  xdot, ydot, thetadot, force, phi, delta, v, vdot, phidot, . . . 
wresid, wdot,  wksi, w, If orce, sf orce, um,  undot)  ; 

(mqq,mqn,  fq,  fc^]=lrgcof  (wksl,  rlocal,rTphi,um,w,  1  force, wresid, . . . 
wdot,  umdot,  grawec, ma,  Im,  n,  mcjq,  mqn,  f q,  mqql,  tarp?,  mcyil ,  . . . 
tarpS,  fql,  teiip9,  tenplO,  Ijnn,  cl, c2,  f  1,  f 2,  f  3,  f  4, betal, beta2,  f q2)  ; 

if  flexlble=l, 

[rm,  gn,kn,  fnn]=smlcof  (nphl,w,mqn,inqq,  fq,wdot,  wresid, rTphl2,  stif,  ... 
grawec,  sforce,ma,lm,  n,  Imn,  fnn,gn,kn,mn,ninnl,knl,kn2,  fnl,  fn2,  fn33,  . . . 
innn,  rlocal,cl,c2,  f  1,  f2,  f3,  f 4,betal,beta2) ; 

(u,udot,u2dot]=lntstnl  (a0,al,a2,a3,a4,a5,a6,a7,u,udot,u2dot,  . . . 
tenpl,  taTp2,  tenp3,  terTp4,  teirpS,  itn,  gn,  f  nn,  kn)  ; 

end 


Iq,qdot,q2dot]  =  intlrg(h,aO,  a3,  a6,a7,mqn,mcjq,  fq,q,qclot,q2dot,  . . . 
taTp6,  u2dot) ; 

(X,  y, theta, xdot, ydot, thetadot, V, phi, vdot, phidot,  x2dot,  . . . 
y2dot, theta2dot, v2dot, phl2dot ] =chvar (q,  qdot,  q2dot,  u,  udot,  u2dot ) ; 

if  rlgld=l, 

(delta, dtheta,tenp20]=rigidcontrol  (mqq,  fq2,  force,  theta, roll, ml2,m21,m22,  . . . 
f  11,  f 22,  hi ,  h2,  bl,  b2.  A,  F,  B,  dtheta2dot,  thetadot,  dthetadot,  . . . 
dtheta,  kv,  kp,  tenpl4,  time,  slope,  tenp20) ; 

[tplot,xplot,yplot,thetaplot,xdotplot, count, deltaplot,dthetaplot]  »=gplot  (. . . 
count,  time, X,  y,  theta, xdot, tplot, xplot, yplot, thetaplot, xdotplot, thetadeg, . . . 
delta,  deltaplot,  dtheta, dthetadeg, dtlietaplot , deltadeg) ; 

(ydotplot,  thetadotplot,vplod,phlplot,vdotplot,plildotplot,count]=gplot2  (. . . 
count,  ydot,  thetadot,  v,  phi,  vdot, phidot,  ydotplot,  thetadotplot,  vplot ,  phlplot, . . 
vitbtplot,phidotplot) ; 

count l=count 1+ 1 

If  count  1=0,  . 

foroe=0; 
end 

If  count  1=0, 
tl=X; 
t2=y; 
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t3==theta; 

t4=xdot; 

t5=ydot; 

t6=thetadot; 

t7=v; 

t8=phl; 

t9=vdot; 

tlO=phicbt; 

tll=x2ciot; 

tl2=Y2dot; 

tl3=theta2cbt; 

tl4=v2(±)t; 

tl5=phi2cbt; 

tl6=delta; 

clear  X  y  theta  xdot  yckjt  thetadot  v  phi  vdot  phidot  delta 

clear  fnl  fn2  fn33  fnn  fq  fql  gn  gnl  kn  knl  kn2  1 force  mn  mnn 

clear  mnnl  mcyi  mqnl  mqq  mqfql  q  qdot  q2dot  sforce  teipl  tenp2 

clear  teiTp3  tenp4  tenpS  teirp6  teiip?  terp8  tGJFp9  u  udot  u2dot 

clear  um  umdot  w  wdot  wksi  wresid  couiitl  tenplO 

clear  x2dot  y2dot  theta2dot  v2dot  phi2dot 

clear  mil  ml2  m21  fl  hi  bl  tenpll  tefTpl2  tenplB  teirpfirst  tenplast 

X=tl; 

y=t2;        '" 

theta=t3; 
'  xdot=t4; 

ydot=t5; 

thetadot=t6; 

v-t7; 

pJii==t8; 

vcbt=t9; 

phidot=tlO; 

x2dot=tll; 

iy2dot=tl2; 

theta2dot=tl3; 

v2dot=tl4; 
.  phi2dot=tl5; 

delta=tl6; 

{q, qdot, q2dot,tenp6,  wresid,  wdot,  wksi,  w,  '. . . 

1  force,  sforce, LBii,undot,  fnn,gn,kn,nm,countl,mll,ml2,n\21,  fl,hl,  . . . 
bl]=constants2a(X,y, theta, xdot, ydot, thetadot, x2dot,y2dot, . . . 
theta2dot) ; 

tu,  udot, u2dot,niinl,  knl,  kn2,  fnl,  fn2,  fn33,nTin,gnl,mqq,mqn,  . . . 
fq,iT*qfll,tenpl,tenp2,tenp3,  tenp4,tenp5,  tenp7,mqnl,tQip8,  fql,taTp9,  . . 
tenplOJ   -  con3tants2b{v,phi, vdot, phidot, v2dot, phi2dot) ; 
end 


end 
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[plotl,  plot2,plot3,plot4, plots,  plot6)=allplotl(tplot,xplot,yplot, . . . 
thetaplot,xckDtplot,yclotplot,tlietadotplot,dthetaplot) 

(Plot7,plot8,plot9,plotl0,plotll]=allplot2{tplot,vplot,phlplot,... 
vdotplot,plilciotplot,cteltaplot) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  lEVEL  2  % 

%       This  function  determines  constant  values  and  initial  values  of  the  % 

%  variables  and  matrices.  % 

%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  (X,y,  theta,xdot,ydot,  thetadot,lin,betal,beta2,E,  Izz,  force,  . . 
phi,  delta, v, iota, epsilon, f time, h, grav,ma, vdot, phidot, n, iim, flexible, 
count, nphl,rrphi2,stif,grawec,rlocal,x2dot,y2dot,theta2dot,v2dot,  . . . 
phi2dot]   =  constantsl 

area=0. 0113097; 

li(t=4; 

b3tal=4.712398/lin; 

beta2=7.853981/ljn; 

count=0; 

delta=0; 

E-2,0e+ll; 

epsilon=0.25; 

flexible=l; 

forcb=30000; 

ftime=0.4; 

grav=-9.8066; 

h=0.001; 

lota=0.5; 

izz=0. 000010178; 

mu=7861..05; 

n=10; 

lmn=lin/n; 

nB=-area*nij; 

p}ii=0; 

ptiidot=0; 

theta=pi/4; 

thetadDt=0; 

X=0; 

xdot=0; 

v=0; 

vdot=0; 

ydot=0; 

y=0; 

rphl=zeros (3,2) ;* 

rTplii2=zeros  (3, 2) ; 

stif=zer63 (3) ; 
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stif (3,3)=E*izz; 

gravvec=zeros (3, 1) ; 

grawec  (3)=grav; 

rlocal=zero3 (3, 1) ; 

x2dot=0; 

y2dot=0; 

theta2dot=0; 

v2dot=0; 

phi2dbt=0; 


%l 

*  I 

%   This  function  determines  constant  values  and  Initial  values  of  the      % 

%  variables  and  matrices. 

%  * 


function  liTi22,  f22,h2,b2,  A,F,B,dtheta2dot,dthetadot,dtheta,kv,  . . . 
kp,  rigid,  tertpl4,  slope,  tenp20]  =constantsla 

ni22=0; 

f22=0; 

h2=0;  r 

b2=0; 

F=0; 

B=0; 

dtheta2dot=0; 

dthetadot=0; 

dtheta=0.7854;%45  Degree 

kv=6; 

kp=9? 

rigid=l;- 

terpl4=zeros  (2) ; 

slope=-2.6180; 

tefp20=0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%   This  function  determines  constant  values  and  initial  values  of  the      % 
%  variables  and  matrices.  % 

%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  (q,  qdot ,  c^dot ,  teitp6,  wresid,  wdot,  wksi ,  w, . . . 

Iforce,  sforce,um,umdot,fnn,gn,  kn,rnn,countl,mll,ml2,m21,  fll,hl, 

bl )  =constants2a  (X,  y,  theta,  xdot,  ydot,  thetadot,  x2dot ,  y2dot,  theta2dot) 
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q^zero3(3, 1) ; 
qdot=zeros  (3, 1) ; 
q2dot=zero3 (3, 1) ; 
teirp6=zero3  (3, 1)  ; 
wresicNzeros  (3) ; 
wctat=zero3  (3) ; 
wksi=zero3  (3, 9) ; 
w=zero3 (3) ; 
lforce=zeros (3, 1) ; 
sforx:e=zero3  {2, 1) ; 
um=zero3(2,l) ; 
uirdot=zero3  (2, 1) ; 
fiin=zero3(2,l) ; 
gn=zero3 (2) ; 
kn^zero3(2) ; 
mn^zeros (2) ; 
count 1=0; 
q(l)=«; 
q(2)=y; 
q(3)=tlieta; 
qdot  (l)=xdot; 
q:tot(2)=ydot; 
q±)t  (3)  =thetadot; 
q2dot(l)=x2cbt;  ,- 

q2dot(2)=y2ckDt; 
q2dot  (3)=theta2dot,• 
mll=ze^os  (2) ; 
ml2=zero3(2,l) ; 
m21=zeros(l,2); 
fll=zeros(2,l); 
hl=zeros(2,l); 
bl=zeros (2,1); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%       This  function  deterndnes  consteint  values  and  initial  values  of  the  % 

%  variables  and  matrices.  % 

%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  lu,ucibt,u2dot,nTnnl,knl,kn2,fnl,fn2,fn33,nnn,gnl,mqq,mqn, 
fq,mcjql,tenpl,terrp2,terp3,  tenp4,  teirpS,  tenp7,mqnl,tenp8,  fql,tGnp9, 
terplO,fq2J  =  constants2b(v,phi,vdot,phicbt,v2cbt,phi2dot) 

u=zeros(2, 1); 
udot=zeros (2, 1) ; 
u2dot=zeros (2, 1) ; 
fnnnl=zeros  (2) ; 
knl=zero3(2); 
kn2=zeros  (2) ; 
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fnl=zeros(2,l); 

fn2=zeros(2,l); 

fn33=zero3(2,l); 

nBin=zero3(2); 

gnl=^ero3  (2) ; 

nKfj=zero3  (3)  ; 

mcji=zeros(3,2); 

fq=2eros(3,l); 

mqql=zero3  (3)  / 

terpl=zeros(2,l); 

tenp2=zero3(2,l); 

taip3=zeros(2,l); 

te!ip4=zeros  (2, 1)  ; 

tetp5=zero3(2,l); 

tarp7=zero3  (3) ; 

mqiil=zero3(3,2); 

teip8=zeros  (3) ; 

fql=zeros(3,l); 

fq2=zeros(3,l); 

tarp9=zeros  (3)  ; 

taipl(>=zeros  (3) ; 

u(l)=v; 

u(2)=plii; 

udot(l)=vdot;  /- 

udot{2)=p}ilcbt; 

u2ciot(l)=v2dot; 

u2ctot(2)=phi2cbt; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%   This  function  determines  the  Initial  values  of  the  output  matrices.      % 
%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  (tplot, xplot, yplot, . . . 

thetaplot,  xdotplot, ydotplot, thetadotplot, vplot, phiplot, vdotplot, 
phictotplot,  thetacieg,deltaplot,dthetaplot,dthetadeg,deltadeg]  =  . 
constants3 (f time, h) 


tplot=2eros (1, (ftlme/h) +1) ; 
xplot=zeros(l, (ftime/h)+l); 
yplot=zeros (1, (ftlme/h) +1) ; 
thetaplot=zero3 (1,  (ftlme/h) +1) ; 
xdotplot=zeros (1, (ftlme/h) +1) ; 
ydotplot=zero3 (1,  (ftlme/h) +1) ; 
thetadotplot=zeros (1, (ftlme/h) +1) ; 
vplot =zeros(l, (ftlme/h) +1) ; 
phlplot=zero3 (1, (ftlme/h) +1); 
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vdot^lot=zeros(l,  (ftljne/h)+l)  ; 

phidotplot=zero3  (1,  (ftJjne/h)  +1)  ; 

thetadeg=0; 

c]eltaplot=zeros  (1,  (ftlnie/h)  +1)  ; 

dthetaplot=zeros(l,  (ftljne/h)+l) ; 

cll;hetadeg=0; 

deltadeg=0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%   Formulating  the  coefficients  used  in  the  sequential  integration  method.  % 
%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  {aO,  al,a2,a3,  a4,a5,  a6,a7]  =  ccoef  (^silon,  iota,  h) 

aO=l/ (epsilon* (h*h) ) ; 

al=iota/ (epsilon*h) ; 

a2=l/ (epsilon*h) ; 

a3=0.5/epsilon-l; 

a'1=iota/epsilon-l; 

a5=h/2*(a4-l); 

a6=h* (1-iota) ; 

a7=iota*h; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
.%  % 

%   Formulating  the  coefficients  used  in  the  mode  shape  function  matrix.     % 
%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%^%%%%%%%%%%%%%%% 


function  lcl,c2,  fl,  f2,  f3,  f 4 J  =TTphicons  (lm,betal,beta2) 

cl=  (sin  (betal*lm)  -sinh  (betal*lm) )  /  (-cos  (betal*lm)  +cosh  (betal*lm) )  ; 

c2=  (sin  (beta2*ljn)  -sinh  (beta2*lm) )  /  (-cos  (beta2*lm)  +cosh  (beta2*lm) )  ; 

e=2*beta2-2* (c2/cl) *betal; 

f 1=1/ (2*cl) + (c2*betal) / ( (cl^2) *e) ; 

f2=-(c2/(cl*e)); 

f 3=- (betal/ (cl*e) ) ; 

f'l=l/e; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%   This  function  calculates  the  matrices  which  changes  with  time.         % 
%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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function  (w,wdot,wksl,wresid,  1  force,  sforce,um,umdlot]  =  . . . 
winatderlwarlant  (theta,X,  y,  x(±)t,ydot,thetaclot,  force, phi,  delta,  v,  — 
vdot,phldot,wrBsld,wdot,wksl,w,  1  force,  sforce,um,umdot) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%   Conputatlon  of  the  transfonnation  matrix  w  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

w(l,l)=l; 

w(2,l)=X; 

w(2,2)=cos(theta); 

w(2,3)=-3ln(theta); 

w(3,l)=y; 

w(3,2)=sin(theta); 

w  (3, 3)  =^03  (theta)  ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%   Cotrputation  of  the  time  derivative  of  the  transformation  matrix  WDCT    % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

wdot(2,l)=xcbt; 

wdot (2,2)=-thetadot*sin(theta) ; 

wdot (2,3)  =-thetadot*cos (theta) ; 

wdot  (3,1)  =ycbt; 

wdot (3,2)=thetadot*cos (theta) ; 

wdot'(3, 3)  =-thetadot*sln  (theta) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%   Cooputation  of  the  derivative  of  the  transformation  matrix  with  respect  % 

%  to  X,  y,  theta:  WKSI  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

wksi(2,l)=l; 

wksi(3,4)=l; 

wksi(2,8)=-sln(theta); 

wksl(3, 8)=co3(theta);  . 

wksi (2, 9)=-co3 (theta) ; 

wksi (3, 9) =-sln (theta) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   Conpilitation  of  the  the  residual  acceleration  matrix  WRESID  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

wresid (2, 2) =- (thetadot"2) *cos (theta) ; 

wresid  (2, 3)  =  (thetadot'"2)  *sin  (theta)  ; 

wresid  (3, 2)-- (thetadot"2)  *sin  (theta) ; 

wresid  (3, 3)  =-  (thetadot"2)  *co3  (theta)  ; 


%   Corrputation  of  the  large  motion  force  matrix 


1  force  (l)=force*  (cos  (theta)  -delta*3ln  (theta)  -phi*sin  (theta) )  ; 
1  force  (2)  -force*  (sin  (theta)  +delta*co3  (theta)  +phl*co3  (theta) )  ; 
Iforce  (3)  =-force*v; 
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%   Conputatlon  of  the  small  motion  force  matrix 


sforce  (1)  =force*  (phi+delta)  ; 
sforce  (2)  =v*foroe; 


%   Carputation  of  the  small  motion  position  vector 


uin(H=v; 
um(2)=phl; 


%   Conputatlon  of  the  small  motion  velocity  vector 


urndot(l)=vcbt; 
uindot(2)=phidbt; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%   This  function  calculates  values  for  the  large  motion  coefficient  matrices% 
%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


fvjnction  [mqq,mqn,  fq,  fq2]=lrgcof  (wksl,  rlocal,nphi,um,w,  Iforce,  . . . 
.wresid,  wdot,umdot,grawec,ma,  Im,  n,mqq,mqn,  fq,mqql,  tenp7,mqnl,  . . . 
tetp8,  fql,  teiTp9,  tenplO,  Imn,  cl,  c2,  f  1,  £2,  f 3,  f  4, betal,beta2,  fq2) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"%%%%%%%%%%%%%%%% 
%   Creating  the  coefficient  matrix  for  M3Q  % 

for  i=l:3, 

terpl=wksi ( : , 3* (i-1) +1 :3*i) ; 
for  j=l:3, 
teTp2=wksi  (:,  3*0-1) +1:3*  j); 
teTp7=tenpl '  *teftp2; 

mqql  (1,  j)=quaclmat  Cast  ring  11',  1/  1/n,  Imn,ma,cl,c2,betal,beta2,  fl,  f2,l3,  . . . 
f4,nphi,0,0,rlocal,tenpl,terp2,un,0,0,0,0,tenp7,0,0,  . . . 
6,0,0); 
end 
end 
nKf:f=na*  OnqRl) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%^%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%      Creating  the  coefficient  matrix  for  MJ^l  % 

%%%%%%%%%%%%%%%%U%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for  1^1:3, 

teitp=wk3l  ( :,  3*  (i-1) +1:3*1)  ; 
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tenp8=taTp  *  *w; 

mqnl  (1,  : )  =quadmat  ( 'astring22  ',1,2,  n,  Imn, ma, cl,  c2, betal, beta2,  f  1,  f 2,  f 3,  . . . 

f  4,nphl,  w,  0,  riocal,  0, 0,  um,  tenp,  0, 0, 0, 0,  tenp8, 0, 0, 0, 0)  ; 
end 
nKy»=nia*  (mqnl);  , 

%       Creating  the  coefficient  mayrix  for  FQ  ^ 


%Cotiputation  of  the  FQ  coefficient 
for  i=l:3, 

taip^ksi  ( : ,  3*  (i-1) +1 :  3*i) ; 

teofD^^tarp'  *wresid; 

fj^ij^^t  CastringSBM,  l,n,lWn,ma,cl,c2,betal,beta2,  f  1,  f2,  f3,  f4,  . . . 

irphi,  0,  wresid,  rlocal,  0, 0,  um,  tenp,  wdot,  umdot,  grawec,  0, 0,  taipS,  . .  - 

terplO,0,0); 
end 

fq2=tTB*fql; 
fq=fc^+lforce; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%   lliis  function 'calculates  values  for  small  notion  coefficient  matrices.   "s 
%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%l%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-6 


.  function  (mn,gn,kn, fnn]=smlcof (nphi,w,mqn,mqq, fq, wdot, wresid,  nphi2, .. . 
stif,  grawec,  s force, ma,  Im,  n,  linn,  fnn,gn,  kn,mn,mnnl,knl,  . . . 
kn2,  fnl,  fn2, fn33,mnn, rlocal,cl,c2, f 1, f2, f 3, f 4,betal,beta2) 

%   Creating  the  coefficient  matrix  W  °, 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%l 
mnn=<5uacknat  ('astring44 ',2,2, n,lmn, ma, cl,c2, betal, beta2,  fl,  f2,f3,  f4, .'. . 

nf^i,  w,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)  ; 
mnn  l=?na  *nnn  ; 
nin=ninnl-niqn  •  *inv  (mqq)  *mqn; 

%   Creating  Uie  coefficient  matrix  GN  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
gnl=^quadknat  (•astring55',2,2,n,  Imn, ma, cl,c2, betal, beta2,  fl,  f2,  f3,  f4, . . . 

nphi,  w,  0, 0, 0, 0, 0, 0,  wdot,  0, 0, 0, 0, 0, 0, 0, 0) ; 
gn=ina*gnl; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   Creating  the  coefficient  matrix  KN  % 

knl=quacknat  ( 'astring66' , 2, 2, n,  ljnii,ma, cl, c2, betal, beta2,  f  1,  f 2,  f 3,  f  4, . . . 


80 


iiphi,  w,  wresld,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)  ; 
kii2=<juadkiet  ('astrlng77',2,2,n,  Iirai,ma,cl,c2,betal,beta2,  fl,  f2,  f3,  f  4,  . . . 

0, 0, 0, 0, 0, 0, 0,  0, 0, 0, 0, 0, 0, 0, 0,  rTphi2,  stif )  ; 
kji=Tna*knl+kn2; 

%%%%%%%%%%%%%%%%%%%%%%%%%%'i%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%''6^ 

%   Creating  the  coefficient  matrix  EN  % 

%l%%%%%%%%%%%%%%%%t%%%%%%%%%%%%%%%%%%%%%%%%%'l%%%%'l%%%%%%%%%%%%%%%%1il%'i%1s%%%%'fel; 
fnl=quadknat  (*astring88',2,  l,n,  Imn,ma,cl,c2,betal,beta2,  fl,  f2,  f3,  f4,nphi,w, . . . 

wresid,  rlocal,  0, 0, 0, 0, 0, 0,  grawec,  0, 0, 0, 0, 0, 0)  ; 
fn33=nia*{fnl)  +  sforce; 
f  nn=fn33-fnqn'  *inv  (mqq)  *fq; 


%  * 

%   This  function  lnpllcltly  integrates  the  small  motion  equation.  % 

%  * 

^ 


function   {u, udot, u2dot ] =lntsml  (aO, al, a2, a3,  a4,  a5,  a6,  a7, u,  udot, . . . 
u2dot,  tenpl,  tenp2,  tenp3,  tejTp4,  teiTp5,mn,gn,  fnn,  kn) 

teiipl=a0*u+a2  *jxtot+a3  *u2dot  ; 
te^p2=al  *u+a4  *udot+a5*u2dot ; 
t  efTp3=nin  *  teirpl ; 
tenp4==gn*  tenp2 ; 
tenpl=f  nn+teiTp3+teitp4  ; 
teiTp2==u; 

tenp5=kn+a0*mn+al*gn; 
terTpl=tenp5\  tenpl ; 
u^terrpl; 
t9ip3=u2dot; 

u2dot=aO* (u-tenp2) -a2*udot-a3*tQip3; 
udot=udot+a6  *  taTp3+a7  *u2dot ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1;%%%%%%%%%%%l%%%%%%'67 
%  ^. 

%       This  function  e>^llcltly  Integrates  tlie  large  motion  equation.  " 

%  ? 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%l%%%%%%%%%%%^ 


function   (q,qdot,q2dot]=lntlrg(h,aO,  a3,  a6,a7,mqii,niq:5,  fq,q,qdot, ... 
q2cbt,  tenp6,  u2clot) 

q=q+h*qdot+a3*q2dot/a0; 
qdot=qdot +a 6 *q2dot ; 
tenp6=tenp6+mqn*u2dot; 
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fcj>=fq-tQtp6; 
q2dot=inv  (inqq)  *fqn; 
qclot=qdot4q2dbt*a7; 
q=qH32cbt/.aO; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%i%%%%%%%%%'i".''u:i 
%  % 

%   This  function  updates  the  tine  dependent  values.  % 

%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'i-.''6 


function  (X, y,  theta,xdot,ydot,  thetadot,v,phl,vdot,plildot,x2dot,  .. . 
y2dot,  theta2dot,  v2dot,phl2dot)  =chvar  (q,  qdot,  q2dot,  u,  udot,  u2dot) 

X=q(l); 

y=q(2); 

theta=q<3); 

xdot=qdot  (1) ; 

ydot=qdbt(2); 

thetadot=qdot  (3)  ; 

x2dot=q2dot(l); 

y2dot=q2dot(2);/-,~ 

tlieta2dot=cj2dot  (3)  ; 

v-u(l); 

ptil=u(2); 

vdot=udot(l); 

phidot=udot(2); 

v2dot=u2dot(l); 

phl2dot=y2dot (2) ; 

%     • 

%   This  function  controls  rigid  or  flexible  missile  with  rigid-body       I 
%  controller.  The  control  angle (delta)  is  limited  to  10  degrees.  " 

%  ".> 


function  (delta,  dtheta,tenp20 J  =rlgldcontrol  (mqq,  fq2,  force,  tlieta,mll,ml2,  . 
m21,m22,  fll,f22,hl,h2,bl,b2,A,  F,B,dtlieta2dot,  thetadot,dthetadot,  ... 
dtheta,  kv,  kp,  tenpl4,  time,  slope,  tenp20) 

if   (time  >=  0)   S    (time  <  0.1), 

dtJieta2dot=0; 

dthetadot=0; 

dtlieta=pi/4; 

tG!Tp20=dtheta; 
elseif  time  >=  0.1, 

dtl»eta2dot=0; 
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dtl  le  tadot=s  Icpe ; 
dtheta=tenp20+slope*  (time-O.l)  ; 
erri 

mll(l,l)=iiiT!(l,l); 

mll(l,2)=mna,2); 

ml  1(2,1)  =Tiiqq (2,1)  ; 

mll(2,2)=inqq(2,2); 

ml2(l)=?Tiqq(l,3); 

ml2(2)=nicn(2,3); 

m21(l)=iTiqq(3,l); 

tn21(2)=inqq(3,2); 

m22=mqq(3,3); 


fll(l)=fq2(l); 
fll(2)=fq2(2); 
f 22=fq2 (3) ; 

hi  (l)=force*cos(theta);  ~^ 

hi  (2) =f orce*sin (theta)  ; 

bl  (l)=-force*sin.(tlieta) ; 
bl  (2)=force*cos  (theta)  ; 
taTpl4=inv(mll) ; 
A=n\22 -fn2 1  *  tenpl  4  *ml  2  ; 
F=f22-m21*tenpl4* (f 11+hl) ; 
B=m21*teTpl4*bl; 

delta=-lnv  (B)  *  (A*  (dtheta2dot-)cv*  (thetadot-dthetadot)  -kp*  (theta-. . . 
dtheta) ) -F) ; 

if  delta  >=  0.174532  ,  %10  degree 

delta=0. 174532; 
Old 
if  delta  <=  -0.174532, 

delta=-0. 174532; 
elseif  (delta  <  0.174532)  |  (delta  >  -0.174532), 

delta=delta; 
efid 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%   This  function  calculates  the  elements  of  matrices  for  output  % 

%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  Itplot,xplot,yplot,tJietaplot,xdotplot,count,deltaplot, 
dthetaplot J  =gplot  (count,  time, X, y,  theta, xdot,  tplot ,  xplot,  yplot, . 
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thetaplot,xdotplot,  thetadeg, delta, deltaplot,dtheta,dthetadeg, . . . 
dthetaplot,  deltadeg) 

count=count+l 

thetadeg= (180*theta) /pi; 

dthetadeg=(180*dtheta) /pi; 

deltadeg= (180*delta) /pi; 

tplot  (count)  =tljTie; 

xplot (count) =X; 

yplot( count )=y; 

thetaplot (count )=t he tadeg; 

xdotplot (count)  =^dot; 

deltaplot (count) =deltadeg; 

dth^aplot  (count)  =dthetadeg; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%       This  function  calculates  the  elemsnts  of  matrices  for  output  % 

%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function    [ydotplot,thetadotplot,vplot,phiplot,vdotplot,phidotplot,count)=. . . 
gplot2  (count,  ycfc)t,thetadot,v,plii,vdot,plildot,ydotplot,thetadotplot,  . . . 
vplot,  phlplot,  vdotplot,  ptiidotplot) 

ydotplot (count) =ydot; 
thetadotplot (count) =thetadot; 
vplot (count) =v; 
phiplot (count)  =phi; 
vdotplot (count)  =vcbt ; 
phidotplot (count)  =phidot; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  % 

%   This  function  plots  the  output  graplis.  % 

%  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  Cplotl,plot2,plot3,plot4, plots, plot6]=allplotl(tplot,}<plot,yplot, 
thetaplot/ xdotplot, ydotplot, thetadotplot, dthetaplot) 

plotl=plot (tplot,xplot) , pause 

plot2=plot (tplot, yplot) , pause 

plot3=plot (tplot, thetaplot, '  — ' , tplot, dthetaplot, '  + ' ) ,  pause 

plot4=plot  (tplot,  xdotplot) , pause 

plot5=plot  (tplot,  ydotplot) , pause 

plot6=plot  (tplot,  thetadotplot) ,  pause 
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%%%%%"u%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%%%%%**%%%%%%% 

%  % 

%   This  function  plots  the  output  graphs.  % 

%  % 


function  fplot7,plot8,plot9,plotl0,plotll]=allplot2 (tplot,  vplot^phiplot, . . . 
vdotplot ,  frfiidotplot ,  deltaplot ) 

plot7=plot  (tplot,vplot) , pause 
plotST^lot  (tplot,phlplot) ,  pause 
plot9-plot  (tplot,vcbtplot) ,  pause 
plotlO-plot  (tplot,phidotplot) ,  pause 
plotl]=plot (tplot, deltaplot) , pause 

%%%%%"6%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  lEVEL  3  % 

%  TiAs  function  will  integrate  the  coefficient  matrices  of  afnall  and  large  % 
%  motion  over  the  length  of  the  missile  using  Sinpscxi's  integration  method.  % 
%  ^  % 

%%%%%"o%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  aa  =  quadrat  (F,   1,  k,  n,   int,  ma,cl,c2,betal,beta2, . . . 
fl,f2,f3,  f4,nphi,w,wresid,  rlocal,teirpl,tatp2,um,tenrp,vdot,umdot, . . . 
grawec,  tenp7,  teiqp8,  tertp9,  tenplO,  rnphi2,  stif ) 


tenpll=zeros (l,k*n) ;  j 

teiTpl2=zeros  (1,  k)  ; 

tarpl3=2eros  (1,  k) ; 

tenpfirst=zeros  (l,k) ; 

tenplast=zeros (l,k)  ; 

aa=zeros (1, k) ; 


x=0; 

for  i=l:n, 

x=x+int; 

terpll  ( : , k*  (i-1)  +1 : i*k)  =feval  (F,x, cl, c2,betal,beta2,  fl,  f 2,  f3, f 4, . . 

nphi,w,wresid,rlocal,tenpl,tenp2,um,tenp,rloccd,wdot,uTdot, . .. 

grawec,  te«p7,terp8,tenp9,ternpl0,nphi2,  stif) ; 

end 

for  m=  l:k, 

for  1  =  2:2: (n-1), 
te[Tpl2(:,m)  =  tenpl2(:,m)   +taipll(:,  k*(i-l)^; 

end 
axL 
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for  p=  l:k, 

for  j  =  3:2: (n-1), 
tefipl3(:,p)   =  terpl3(:,p)    +tenpll(:,   k*(j-l)+p); 

end 

tcjipf irst  (: ,  1  :k)  =tenpll  ( : ,  1  :k)  ; 
tenplast  (,: ,  1  :k)  =tenpll  ( : , n*k-  (k-1)  :n*k)  ; 

aa=(lnt/3)  *  (tenpfirst+4*tenpl3+2*tenpl2+tenplast)  ; 


-o  «%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


9-0. ' 

'^'^WUUWU««««eoooooooaag0U0OO00e00O0OOOQD0D0000000DDQD00D0000D0000D0000000D0004 

%  % 

%       This  function  creates  the  product  of  matrices  for  integration  of  M3Q  % 

%  coefficient  matrix.  % 

%  % 


function  yll=astrlngll  (x,  cl,c2,betal,beta2,  fl,f2,  f3,  f4,rTphl,w, . . . 
wiGsid,  rlocal,terTpl,tefTp2,imi,  tenp,  rlocal,wdot,umci3t,gravvec,  teirp?, , . . 
tatpa,  teiTp9,  tarpl0,nphi2,  stif ) 

rrj'hi  (3,  l)=fl*  (cl*  (cos  (betal*x)  +cosh  (betal*x) )  +sin  (betal*x)  +. . . 
s  inh (betal*x) ) +f3*  {c2* (cos (beta2*x) +cosh (beta2*x) ) +. . . 
sin (beta2*x) tsinh (beta2*x) ) ; 

nphi (3, 2)=f2* (cl* (cos (betal*x) +cosh (betal*x) ) +sin (betal*x) +. . . 
si  nh (betal*x) ) +f 4* (c2* (cos (beta2*x) +cosh (beta2*x) )  +  . . . 
sin (beta2*x) +slnh (beta2*x) ) ; 

r]ocal(l)=l; 
rlocal(2)=x; 

yl]=rlocal'*teiip7*rlocal+rlocal'*tenp7*rTphl*unH-um'  *nplii'*teirp7*. . . 
nphi*um+um' *nphl  •  *tenp7*rlocal; 

?;"o''6%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%       This  function  creates  the  product  of  matrices  for  integration  of  M3^  % 

%  coefficient  matrix.  % 

%  % 

%°6''6%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%*%%%%%*%%*%*%% 


function  y22=astrlng22(x,cl,c2,betal,beta2,fl,f2,f3, f4,nphi,w, ... 
wresid,  rlocal,  tenpl,  taTp2,  urn,  tenp,  rlocal,  wdot,  umdot, grawec,  tenp7,  . . . 
tenpB,  tenp9,  tenplO,  nphl2,  stif) 

nplii  (3,  l)-f  1*  (cl*  (cos  (betal*x)  +cosh  (betal*x) )  +sln  (betal*x)  +. . . 
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si  nil  (betal*x) )  +f  3*  (c2*  (cos  (beta2*x)  +cosh  (beta2*x) )  +. . . 
sin (beta2*x) fsinh (beta2*x) ) ; 

rtplil  (3,2)=f2*(cl*(cos(betal*x)+cosh{betal*x))+sln(betal*x)  +  ... 
sinh (betal*x) ) +f4* (c2* (cos (beta2*x) +cosh (beta2*x) )  +. . . 
sin (beta2*x) +sinh (beta2*x) ) ; 

rlocal(l)=l; 
rlocal(2)=x; 

y22=rlocal  •  *tenp8*nphi+um'  *nphi  •  *tenp8*iTphi; 

%%%!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%       Tills  function  creates  the  product  of  matrices  for  integration  of  FXJ  % 

%  coefficient  matrix.  % 

%  % 

%%''b1,%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  y33=astring33(x,cl,c2,betal,beta2,fl,f2,f3,  f4,nphl,w, 

wrcsid,  rlocal,teirpl,tefTp2,um,  terp,  rlocal,  wdot,umdot,gravvec,  teirp?, . . . 
teiip8,  tefrp9,  tenplO,  nphi2,  stif ) 

nphi (3, 1) =fl* (cl* (cos (betal*x) 4cosh (betal*x) ) +sln (betal*x) +. . . 
sinh (betal*x) ) +f 3* (c2* (cos (beta2*x) +cosh (beta2*x) ) +. . . 
sin (beta2*x) +sinh (beta2*x) ) ; 

nphi (3, 2) =f2* (cl* (cos (betal*x) +cosh (betal*x) ) +sin (betal*x) +. . . 
sinti  (betal*x) )  +f  4*  (c2*  (cos  (beta2*x)  +cosh  (beta2*x) )  + . . .  j 

sin (beta2*x) +3inh (beta2*x) ) ; 

rlocal(l)=1.0; 
rlocal(2)=x; 

y33=-rlocal '  *teitp9*rlocal-rlocal '  *tefTp9*nphi*um-2*r local '  *teiTplO* . . . 
apl  I  i  *unidot-um'  *iTphi '  *tenp)9*rlocal-um'  *rTfihi  *  *tenp9*nphl*um-2  * . . . 
um '  *rTphi '  *tenplO*nplil*umdot+rlocal  *  *t€iTp'  *grawec+grawec '  *tenp* . . . 
tiphi*Lin; 

%%"^%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%       This  function  creates  the  product  of  matrices  for  integration  of  m  % 

%  coefficient  matrix.  % 

%  % 

%%%'!;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function  y44»=astrlng44  (x,cl,c2,betal,beta2,  fl,  f2,  f3,  f4,nphl,w, . .. 
wresld,  rlocal,  tenpl,  tenp2/Um,  tenp,  rlocal,wdot,uiTciot,grawec,  teftp7,  . . . 
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tenp8,  tenp9,  terplO,  rrpJil2,  stl  f ) 

irphi  (3,  l)=f  1*  (cl*  (cos  (betal*x)  ^cosh  (betal*x) )  ^sln  (betal*x)  +  . . . 
si  nh (betal*x) ) +f 3* (c2* (cos (beta2*x) +cosh (beta2*x) ) + . . . 
sin (beta2*x) +sinh (beta2*x) ) ; 

iiphl  (3,2)=f2*  (cl*  (cos  (betal*x)  +cosh  (betal*x) )  +sln  (betal*x)  +  . . . 
sinh (betal*x) ) +f 4* (c2* (cos (beta2*x) +cosh (beta2*x) )  +. . . 
sin (beta2*x) +sinh (beta2*x) ) ; 

v'l'l=Tiphi '  *w'  *w*nphi; 

°c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%   This  function  creates  the  product  of  matrices  for  Integration  of  GN     % 
%  coefficient  matrix.  % 

%  % 

°5^%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  y55=astring55(x,cl,c2,betal,beta2,fl,f2,  f3,  f4,nphi,w, . . . 
tviesid,  rlocal,  tenpl,  tenp2,um,  tenp,  rlocal,wdot,unidot,grawec,  tenp7, . . . 
tenp8,  tenp9,  teiTplO,nphi2,  stif ) 

nphi (3, 1) =f 1* (cl* (cos (betal*x) +cosh (betal*x) ) +sin (betal*x) +. . . 
sinh (betal*x) ) +f3* (c2* (cos (beta2*x) +cosh (beta2*x) ) +. . . 
s  In (beta2*x) +sinh (beta2*x) ) ; 

nphi (3,2)=f2* (cl* (cos (betal*x) +cosh (betal*x) ) +sln (betal*x) +. . , 
Einh(betal*x) ) +f4* (c2*  (cos (beta2*x) +cosh(beta2*x) ) +. ..        ^ 
sin (beta2*x) +sinh (beta2*x) ) ; 

y55=2*nphi  •  *w  •  *wdot*nphi  ; 

%  % 

%   This  function  creates  the  product  of  matrices  for  integration  of  FN     % 

%  coefficient  matrix.  % 

%  % 


function  y66=astrlng66(x,cl,c2,betal,beta2,  fl,f2,  f3,  f4,nphl,w, ... 
wresid,  rlocal,  teiipl ,  tenp2,  urn,  tenp,  rlocal,  wdot, umdot,  grawec,  teiTp7, 
tcjfpB,  terTp9,  tenplO,  nphi2,  stif ) 

nplii  (3,  l)=fl*  (cl*  (cos  (betal*x)  +cosh  (betal*x) )  +sln  (betal*x)  +. . . 
sinh (betal*x) ) +f3* (c2* (cos (beta2*x) +cosh (beta2*x) ) +. . . 
si  n (beta2*x) +slnh (beta2*x) ) ; 
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irphl  (3, 2)  =f2*  (cl*  (cos  (betal*x)  +cosh  (betal*x) )  +sln  (betal*x)  + . . . 
sinh (betal*x) ) +f4* (c2* (cos (beta2*x) +cosh (beta2*x) ) +, . . 
sin (beta2*x) +sinh (beta2*x) ) ; 

y66-Trphl '  *w'  *wresid*rTphl; 

%%"i'6^%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  % 

%   This  function  creates  the  product  of  matrices  for  Integration  of  KN     % 
%  coefficient  matrix.  % 

%  % 

%%?i7.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  y77=astrlng77  (x,  cl,c2,betal,beta2,  fl,  f2,  f3,  fA,Tfphi,w,  . . . 
wrcsid,  rlocal,tenpl,terp2,um,  teitp,  rlocal,wdot,umdot,grawec,  tefrp7, . . . 
teipO,  teTp9,  teiTplO,  nphi2,  stif ) 

rTr'hi2  (3,  l)=fl*  (betal^2)  *  (cl*  (-cos  (betal*x)  +cosh  (betal*x) ) . . . 

-sin  (betal*x)  +sinh  (betal*x) )  +f  3*  (beta2'^2)  *  (c2* . . . 

(-cos (beta2*x) +cosh (beta2*x) ) -sin (beta2*x) +slnh (beta2*x) ) ; 

rphi2  (3, 2)=f2*  (betal'^2)  *  (cl*  (-cos  (betal*x)  +cosh  (betal*x) ) . . . 

-si  n  (betal*x)  +slnh  (betal*x) )  +f  4*  (beta2'-2)  *  (c2* , . . 

(-cos (beta2*x) icosh (beta2*x) ) -sin (beta2*x) +sinh (beta2*x) ) ; 

y77=7rphl2 '  *stif  *iTphi2; 

%%?.'^%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  ./  % 

%   This  function  creates  the  product  of  matrices  for  Integration  of  FN     % 
%  coefficient  matrix.  % 

%  % 

%%7,%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  y88=astrlng88(x,cl,c2,betal,beta2,  fl,f2,  f3,  f4,rrphl,w,  ... 
wresid,  rlocal,  tenpl,  tenp2, urn,  tenp,  rlocal,  wdot, umdot, grawec,  tefTp7, 
teripB,  tenp9,  tenpl  0,  itphl2,  stif) 

rtphi  (3, 1)  =f  1*  (cl*  (cos  (betal*x)  +cosh  (betal*x) )  +sin  (betal*x)  +. . . 
sinh (betal*x) ) +f 3* (c2* (cos (beta2*x) +cosh (beta2*x) )  +. . . 
sin (beta2*x) +sinh (beta2*x) ) ; 

nplii  (3, 2)  =f2*  (cl*  (cos  (betal*x)  +cosh  (betal*x) )  +sin  (betal*x)  +. . . 
si  nh (betal*x) ) +f 4* (c2* (cos (beta2*x) +cosh (beta2*x) ) +. . . 
sin (beta2*x) +sinh (beta2*x) ) ; 

r]ocal(l)=l; 
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rlocal  (2)=x; 

y08=wphl  •  *w'  *gravvec-nphi '  *w'  *wresid*rlocal; 
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